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SUMMARY 


This manual is designed to explain the intent and broad outline 
of the computer program FRAC3D, which is designed for use in analyzing 
stresses in structures (including plates, bars, or blocks) which may con- 
tain part-circular surface cracks or embedded circular cracks. Instruc- 
tions are provided for preparing input, including that for the supporting 
programs LATTICE and MATSOL as well as for FRAC3D, and the course of a 
substantial illustrative calculation is shown with both input and output. 
The formulas underlying the calculations are summarized and related to 
the subroutines in which they are used. Many issues of strategy in using 
this program for analysing stresses around surface cracks are elucidated 
by applications of it that are discussed in the report to which this 
manual is a supplement (NASA CR 159400) . 

INTRODUCTION 


This manual describes the capabilities of the computer program 
FRAC3D which is based on the theory developed in [1]*. The analysis capa- 
bility consists of solving for the three-dimensional stress state in a 
three-dimensional region which may or may not contain a flaw. The types 
of flaws analyzable by the program include completely buried circular flaws 
or part-circular surface flaws. The manual begins with a description of 
the allowable geometric configurations. This section is followed by a brief 
synopsis of the various analytical forms which describe the stresses and 


* Numbers in brackets refer to entries in the list of references. 
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displacements and of the procedure of developing the needed input data and 
their various formats. A representative stress analysis is used as a com- 
prehensive illustration, with discussion of the strategy, the input and 
several forms of output. One appendix summarizes the formulas underlying 
the program, and another appendix lists its component parts, including a 
brief description of each and references to the formulas each employs. 

The program described here is operational on Battelle f s Columbus 
Laboratories 1 CDC Cyber 73 computer system. Typical computer times on this 
system are ten to twenty minutes for a substantial three-dimensional analy- 
sis, with considerable variation possible according to the extent of detail 
employed in the analysis. Many illustrative results are included in the 
report which this manual accompanies [2]. 

RANGE OF APPLICABILITY OF THE PROGRAM 
Geometries Considered 

The geometries analyzable with the computer code FRAC3D range from 
a circular flaw in three-space to a tilted, part-through crack in a right- 
parallelepiped as shown in Figure 1.* Shown there is a six-faced solid, all 
edges meeting at right angles, with a tilted circular surface flaw cutting 
the top surface. Each face is numbered and has associated with it a local 
rectangular coordinate system. The crack is described by the angle w and 
the directed distance labelled r c * The program uses the crack radius a as 
the reference (unit) length. The tilt angle, can be in the range 

- but must be less than -y for the crack to cut the top surface. 

I l l (1) 

The quantity r is the distance between the crack center and the (x , 

(1) (1) C 

y , z ) origin, with sign being positive if the center is outside the 
body or negative if it is inside. 


* The crack also can be omitted, so that the program can be used for 
analyzing stresses in a half-space, quarter-space, slab or block sub- 
jected to a wide variety of boundary conditions. 
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FIGURE 1. COORDINATE SYSTEMS ASSOCIATED WITH A 

RECTANGULAR BLOCK WITH A CIRCULAR CRACK 
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The local rectangular coordinate systems shown in Figure 1 and 
the numbering scheme which identifies them are fixed conventions within the 
program. The cylindrical system associated with the crack arbitrarily has 
been assigned the identification number 7. Global coordinates (x,y,z) are 
taken to coincide with those of the system labelled 1, and the geometry of 
a region to be analyzed is defined relative to this system. Ignoring the 
crack for the moment, the six faces are defined by the planes listed in the 
following table. 


TABLE 1. PLANES DEFINING ANALYZABLE REGIONS 


Face No. Plane Local Coordinate System 


1 

N 

II 

O 

X (1) , 

y (1) , 

z< L > 

2 

x = x A 

* (2) , 

y (2 >. 

z< 2 > 

3 

y = y« 

x (3) . 

(3) 
y > 

z< 3 > 

4 

At 

z = z 

x (4) . 

y< 4 >, 

z< 4 > 

5 

L 

X = X 

x (5) . 

y< 5) . 

z< 5 > 

6 

u 

y = y 

x< 6 >. 

y <6) , 

z (6) 


The user can include or exclude any meaningful combination of the 
seven numbered coordinate systems to describe various regions to be analyzed. 
Thus a half space with a crack can be described by choosing only Face 1 and 
the crack system. Similarly, a finite thickness slab can be described by 
choosing either face pairs 1 and 4, or 2 and 5, or 3 and 6, in addition to 
the crack system. A slab with a crack along an edge might be described by 
Faces 1, 3, and 6 in addition to the crack, and so forth. 

Loads Considered 


The analyses underlying FRAC3D express both stresses and displace- 
ments in the body in terms of the influence of certain constants which de- 
fine loads acting either on the crack or on any of the six faces. If all 
the load constants are known initially, they can be used directly in the 
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"RESULT" mode of operation of the program to find stress and displacement 
components anywhere in the body. If the load constants are not known, but 
instead a definite set of their implied boundary tractions (and/or displace- 
ments) is known, then it is necessary first to use the program in another 
mode ("SOLVE") to find all the load constants. This is the more common way 
in which the program has been used. 

In considering stresses associated with cracks, it is customary to 
take the defining boundary tractions to be the stresses which would be trans- 
mitted across the plane of the crack if no crack were present, but which are 
not transmitted when the crack is there. Such tractions ordinarily have a 
fairly limited range of influence, and they fully determine the stress inten- 
sity factors along the crack front. The stresses and displacements they 
imply can be added to those which would arise in the absence of the crack in 
order to get a full description of the stresses and displacements in the 
body. For some situations, it may be desirable to use the program in a pre- 
liminary calculation to find the stress patterns that would exist in the 
crack location if no crack were there. The program allows evaluation of such 
stresses under a wide variety of loadings. 

ANALYTICAL FORMS AND OPERATIONS APPEARING IN THE PROGRAM 


The intent of this section is only to present typical analytical 
forms of stresses and displacements and to describe in broad terms how in- 
formation flows in a calculation, so that the user can understand the vari- 
ous input and output quantities. Derivations and further discussion are 
given in References 1, 3, and 4. 

Stresses and displacements arise from surface loads applied to the 
crack or to any of the six plane surfaces bounding the body. Load constants 
associated with the crack or with any other single surface are used to 
define those surface loads, but the load constants, for any one surface, 
are not alone in implying tractions on that surface, since loads applied 
anywhere may imply stresses everywhere. A primary goal is to coordinate 
the choice of all the load constants in a way to provide good satisfaction 
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of overall boundary conditions, and to this end, generalized expressions 
are used for stresses from loads on any one surface. 


Expressions Related to Crack Loads 


The expressions for stresses associated with the crack constants 
are shown by Equations (1) through (6) of Appendix A. For example, the 
expression for , the normal stress in the direction perpendicular to the 
plane of the crack, is 


y m J 0 k ^ 0 j[“m,k F m,k,l + B m,k F n/k,2 + Y m,k F mfk,3] 


cos m 0 


[ a m,k F mfk, 1 + P m,k F mjk,2 + Y m,k F mfk,3] sin m6 [ 


Here a ,,3 , , Y , , a ,,3 , , and y . are crack load constants, and p 

m , k. m , X m y k. m y k m , k m y k 

is the shear modulus. Constants of the first three kinds relate to loads 
essentially symmetric around 0=0, while those of the last three kinds 
refer to conditions antisymmetric around 0 = 0 so that they can often be 

i * i 

dropped. The constants a , and a , refer to normal loads on the crack 

m , x m , x 

face, while the other four sets of constants refer to tangential loads on 
the crack face. Thus, for situations in which only normal loads would be 

applied to the crack and in which there is symmetry around 0=0, the only 

» 

relevant load constants are those denoted as a . . The influence functions 

m,k 

for all these constants, as they act to produce stresses , are the 
trigometric factors cos m0 or sinm0 multiplied by one of the functions 
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(p,C) [ I m,m+2k + l + I m,m+2k + lJ , 


F m!k,2 (p,?) 2 I m,m+ 2 k +2 


F ifk,3 (p ^> 


2-v in,m 


il 2 

2 m,m + 2 k 


for k = 0 


for k > 1 


Here, v is Poisson's ratio, p = r/a, £ = y/&, where a is the crack radius, 
and 


(P.C) ■ I e -«l'U n - >s J M (P.E) J N + ^ (Ode . 

0 

Since integrals of this form perforin such an important role in the analysis 
and special means are needed to evaluate them, several subroutines of the 
program are devoted to their evaluation. The mathematics involved in this 
evaluation is discussed in Reference 1. Here it is being observed merely 
that the stress components (and also displacement components) associated 
with crack load constants are sums of series of those constants multiplied 
by influence functions which the program calculates from geometric specifi- 
cations and from the indices identifying the load constants. 

One aspect of treating stresses associated with crack loads is 
that the user must decide which of all the possible crack constants are to 
be used, that is which combination of the indices m and k will be admitted. 
This is not always a simple matter, since the desirable extent of the series 
depends on the amount of detail needed to provide an adequate fitting of the 
stress field throughout the body. The program currently is prepared to em- 
ploy values of m from 0 through 16, and values of k from 0 through 10, but 
not all of the implied 187 combinations of these indices are necessarily 
desirable. In particular, if the part of crack circle penetrating the body 
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is less than a semicircle, if the load is normal, and if there is symmetry 
around the crack plane and around the plane 0=0, then it appears desir- 
able to drop all the odd values of m. The remaining even values of m pro- 
vide series which furnish complete generality for expansions in the half 

TT 7T 

space where - y < 0 <. y > and they avoid the ambiguity which would arise 
from employing also the superfluous odd values of nu Other combinations 
of indices too may be difficult or unnecessary to use, but it has been found 
that good results can be obtained for several types of problems using 9 to 
45 selected crack constants. Reference 2 describes several informative 
cases. Some experience may be helpful with problems of a new type before 
the final choice of crack constants is made. 

A property that makes the crack load constants particularly 
important is that they are the load constants which determine the stress 
intensity factors along the crack front. In particular, for Mode I loading 
(that is only normal loads on the crack) in which there is symmetry around 
0=0, the stress intensity is* 


K 


I 


l 

m=0 


tt 

a cos m 0 , 

m 


ti 

a 

m 




where 


The stress intensity calculated by the program MATSOL, which solves the 
system of equations set up by FRAC3D, is the ratio of this divided by 
the stress intensity factor that would arise in an infinite body with 

a complete circular crack subjected to a unit normal load. Thus, after 
computations presuming any load, to get a dimensional value of the 
computed K^/K^ should be multiplied by 2p^ / a/ir where p^ is a unit stress. 

* This formula presumes the definition 
Kj = / 2 tt (r-a) (r,6,0). 

2 , — 

If a uniform normal load p is used, this implies K T = — p / ira . 

r o 9 Y I °° tt r o 
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Expressions Related to Surface Loads 


The loads associated with any one of the six surfaces are express- 
ed in terms of elemental loads distributed over rectangular areas on that 
surface, each with intensity that varies in the pyramidal form shown in 
Figure 2. Elemental loads such as these, chosen with suitable bases and 
peak values, can be superposed to give efficient and continuous representa- 
tion to widely varying load patterns. To define such a representation, one 
must first select the lattice of rectangles used for bases of the pyramidal 
elements. It should be understood that the lattice need not be regular, 
but enough lines should be used to provide an adequate number of pyramidal 
bases, each being a rectangle subdivided into four rectangles, as shown in 
Figure 2. In a properly drawn lattice, a pyramidal base can be found around 
any point where two lattice lines cross and extend in both ways beyond the 
"pyramid" point, and the peak pyramidal load is the value of the elemental 
load at that point. Three kinds of loads may be applied on any pyramidal 

base on a surface = 0, namely p - - a^, s~t^, and t~x^. The 

z zx y z 

three peak values p , s , and t associated with the particular pyramid 

m m m 

point become part of the overall set of needed load constants which may 
be specified or sought as needed to finish defining the surface load. The 
pyramidal bases are treated as being part of the initially assigned geome- 
try. 

The program incorporates formulas for computing stresses and dis- 
placements that may be caused anywhere by an elemental load with peak value 

p , s , or t acting on a particular base, treating the body as a half space 
m m m 

bounded by the plane of the base. Thus, for example, the contribution to 
a from the elemental loads acting on a particular base is 


a 

z 


zz 

= P m K l 


ZZ 

S m K 2 


ZZ 

+ t KI 
m 3 


where , and are influence functions which the program computes 

from the geometry of the pyramidal base and the location of the point where 
the stress is to be evaluated. To find the influence function, the program 
uses weighted sums of certain elementary functions evaluated nine times — 
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FIGURE 2. TYPICAL FORM OF SURFACE LOAD ELEMENTS 
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each evaluation being based on distances from one of the nine "corner" 

points (lattice points) in the pyramidal base to the point where the stress 

is to be evaluated. Thus, for example, if the stress a is to be evaluated 

z 

at (x,y,z) for a load p applied to the top face, then 

*1* (x,y,x) = - ^[zP + ag arctan ^ j 

/2 2 2 ~ 

where a = x £ - x, $ = - y, p = /a + g + z , and the asterisk denotes 

the weighted summation used in combining evaluations for each of the nine 
pyramidal base corners (x^, y^, 0) . (Details may be found in Appendix A.) 

To perform this kind of operation, the computer must be supplied the coor- 
dinates of nine corner points of each pyramidal base. To facilitate this 
process, an auxiliary program called LATTICE has been prepared which com- 
putes the information required, and identifies it using index numbers 
assigned to each lattice point. The information is punched on cards in an 
arrangement suitable for input to FRAC3D. 

The input to the program LATTICE consists of coordinates of end 
points of each line of the lattice chosen on each of the six plane faces 
that is used. The process of choosing the lattice design is not rigidly 
formulated, but for cases involving a crack, the process can be elucidated 
by preliminary calculations (also from FRAC3D) showing infinite-body 
stresses needing to be freed when the body is finite. Directions for and 
illustrations of these matters are furnished later in this manual. Many 
examples of lattices are furnished in the accompanying report [2]. 

Assembly of Boundary Conditions and Equilibrium Conditions 

The purpose of the program FRAC3D operating in its "SOLVE" mode is 
to construct a set of linear boundary condition equations which can be 
solved by a least-squares process to determine an appropriate set of crack 
and surface load constants. This operation amounts to constructing a matrix 
of influence functions for unit loads of the kinds represented by the various 
unknown load constants acting to produce surface tractions of three kinds at 
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selected surface points. Each column consists of influence functions for a 
particular load constant and each row shows influence functions for all the 
constants acting to produce a particular type of traction at a given bound- 
ary point. In addition, it is possible to append one or several rows to the 
matrix expressing overall equilibrium conditions which should be satisfied 
exactly rather than in the least-square sense. 

The order of appearance of the load constants in this system of 
equations places the surface-load constants for Surface 1 first, then those 
for Surface 2 (if they are used), and so on through Surface 6, and finally, 
the crack-load constants. Among the surface constants, the p, s and t for 
the first pyramid point are placed first, then those for the second pyramid 
point, and so forth. Among the crack constants, those for m = 0 come first, 
being ordered by increasing k, then those for m = 1, and so forth. This 
basic ordering is arranged by the computer through use of the data arranged 
for FRAC3D by the preliminary program LATTICE. The set of load constants 
is frequently reduced later when symmetries are accounted or when certain 
crack constants are withdrawn from consideration, but the ordering of the 
remaining load constants remains the same as in the basic arrangement. 

The ordering of the rows of the matrix depends somewhat on the 
order in which the user assigns boundary conditions. He is required to 
assign those for Face 1 first (if any are to be used for it) then those 
for Face 2, and so on through Face 7 (the crack face). For each assigned 
point, the computer arranges three equations, namely conditions for 

and x^in that order for a point on an outer surface, or a, , T n , 
yz zx > 

and x on the crack, but the order in which the boundary condition points 
are assigned for a given face is chosen by the user. This is an option 
which can facilitate special studies such as of the behavior of influence 
functions or of convergence of final solutions, though changing the order 
of points for a given face should not alter the solution for load constants 
found ultimately by the use of MATSOL. 

The program allows the assignment of boundary conditions at points 
anywhere on any face or on the crack, including the crack front, though 
boundary conditions at the crack tip may be poorly representative for those 
on the front face. Crack functions for points off the crack plane are 
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evaluated In part by a quadrature which loses accuracy significantly for 
points quite near the crack (say for |£| << 0.01 a), but an alternate 
evaluation procedure called automatically when C= 0 is quite accurate and 
thus permits free assignment of boundary conditions on the crack. (Even 
for points near but not on the crack plane another special subroutine pro- 
vides usable values, though a little slowly. Evaluations of crack functions 
for points near the cylindrical axis of the crack coordinate system is 
handled by still another special calculation.) The continuity of the sur- 
face load functions also allows assignment of boundary condition points 
anywhere on the plane surfaces, but this is a freedom that should be used 
judiciously. Separate research into this matter has shown that the proper 
selection of boundary conditions points can be a vital matter. This matter 
is discussed later and in the accompanying report [2]. 

The basic theory of stresses arising from crack loads presumes 
that the tractions effectively acting at any point on one face of the crack 
are equal in magnitude, but opposite in direction to those acting at the 
equivalent point on the other face. Thus, the crack loads should not be 
expected to produce any net force to upset the overall equilibrium of a body 
containing all or part of the circular crack. This implies in turn that 
the surface loads chosen as freeing forces should themselves be in overall 
equilibrium, but this condition is not implied automatically by the surface 
load theory. Therefore, the program has in it a procedure for constructing 
equations which must be satisfied by all the surface load constants being 
used if the body is to be in equilibrium in each of the three principal 
directions and if there is to be no net moment about the three principal 
axes. These equations take into account the size and location of the base 
over which each elemental load acts, and assign to each element a center of 
action so that overall moments can be computed. Thus, as many as six exact 
equations can be appended to the set of boundary conditions, reflecting the 
need of the body to be in equilibrium in each of six senses. Symmetry may 
make some of these equations identically zero, and these should be omitted 
lest they produce a singular matrix. Moreover, a user may simply desire 
some equilibrium conditions to be unenforced, so the user is simply given 
the option of satisfying as many of the six equilibrium conditions as he 
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may desire. For the case of a slab with a particular surface crack situated 
and loaded symmetrically about both the x and y axes, for example, this means 
that ordinarily it is appropriate to enforce one equilibrium equation (for 
the sum of normal loads on the two surfaces), but it may be omitted if the 
user desires. (It might seem that adequate satisfaction of pointwise bound- 
ary conditions on the two surfaces would automatically insure satisfaction 
of this overall equilibrium condition, but in practical calculations that 
expectation is often unfulfilled and direct imposition of the equilibrium 
condition seems beneficial.) 

The call for imposition of overall equilibrium conditions of the 
six possible kinds is inserted on a yes or no basis for each kind of 
equilibrium, as is explained later. 

Allowance for Symmetries 

If the shape of the body and the loading pattern are symmetric 
about any of the planes x = 0, y = 0, or z = z t^ 9 or ^ there is symmetry 
about any of the lines x = y=0, orx=z - z^/2 =0, or y = z - z^/2 = Q s 

then many pairs or even sets of four or eight load constants are predictably 
equal or equal except for sign. In such cases, it is very desirable to 
reduce the number of unknowns to be found by combining the influence of all 
predictably equivalent constants and then solving for only one representa- 
tive member of each set. To this end, options are provided to make the pro- 
gram allow for symmetries of one or more of the six kinds listed above. In 
making such an allowance, the program first sorts through the list of load 
constants to determine which sets of them are predictably equal (except pos- 
sibly for sign), then selects a representative constant from each set and 
combines the influence functions for all members of the set into an overall 
influence function to be multiplied by the representative constant of the 
set. The mechanism for doing this is a chain of operations which accumu- 
lates the necessary information for proper disposition of each load constant 
and combines the influence functions for related load constants as required. 
The symmetries in the roles of the constants are recorded in a list of 
numbers called A(I) in FRAC3D or E(I) in MATSOL. If A(I) = 0, it is implied 
that the I — load constant is one of the chosen representative constants, 
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so it is not to be displaced. For constants which must vanish because of 

symmetry (such as t on a plane of antisymmetry at x = 0), the program puts 

zx 

A(I) = I. If A(I) is shown to be some other positive (or negative) number, 
it is implied that influence functions for the I— load constant should be 
added to (or subtracted from) the influence functions for the load constant 
with position given by the other number. At the end of the reduction process, 
the combined influence functions for those constants having their A(I) = 0 
are retained, but the columns of influence functions for all other load con- 
stants are deleted. 

Symmetry which causes predictable equalities between load constants 
must include symmetry in the boundary conditions applied. Therefore, a 
further convention employed by the program is that boundary conditions will 
be prescribed at only one of any set of symmetrically situated points, it 
being presumed that symmetric boundary conditions are implied at all other 
points of the set. Since application of a boundary condition at any point 
of a symmetric set leads to the same reduced equation (that is the same 
equation connecting the representative load constants), the effect of apply- 
ing conditions at all points of the set can be had in the least squares 
solution by employing the equations for tractions at just one point of each 
set and multiplying them by a weight equal to the square root of the number 
of points in that symmetry set. The program does that also automatically. 
Thus, for a slab with a surface crack and symmetry around the planes x = 0 
and y = 0, boundary conditions applied on the front face at the origin 
(where x = y = z = 0) are assigned weight 1, those elsewhere on an axis 
are assigned weight /T", and those on neither axis are assigned weight 2. 

The symmetry chain manipulates crack constants and their influence 
functions as well as surface load constants and their influence functions. 
Therefore, it offers a convenient approach to editing the solution process 
in further ways that may be desired. For example, it is possible to use it 
to delete the crack load constants for odd m, together with their columns of 
influence functions. This can be accomplished by a strategically situated 
insertion that A(I) = I for all I related to an odd m. 

The means for asserting the kinds of symmetry to be imposed consist 
of numbers the user supplies showing, in general terms, how tractions of a 



16 


given kind in one region are related to those of the same kind in another 
region where symmetry might exist. Details for doing this are provided 
later. It may be added that the basic theory of stress analysis requires 
that along an edge the shear components of stress acting perpendicular to 
the edge on both the faces must be equal. Since this is a universal require 
ment, provisions for satisfying it automatically have been included in the 
program. 


Alternative Modes of Operation of FRAC3D 

When it is desired to find load constants to satisfy a given set 
of boundary conditions on the body, then FRAC3D is used to derive a set of 
boundary conditions by operating it in the "SOLVE" mode as has been 
described. A modification of this mode is had by adding the instruction 
"OUTPUT*', and when this is done, the computer prints out the influence func- 
tions computed for each of the boundary equations. This offers a convenient 
tool for finding, for example, what stresses would exist on a given plane in 
an infinite space having a circular crack subjected to a given kind of load- 
ing (that is loading of a pattern related to specific crack load constants) . 
If the plane is to be a surface of a finite body, then it is possible to 
draw a map of the boundary stresses which are to be freed, except as they 
are modified by interactions between surface and crack loads. Such a map 
is useful in planning a lattice for the surface load functions. 

If the load constants are known, then operation of FRAC3D in its 
"RESULT" mode makes it possible to compute the resultant stress and displace 
ment components anywhere in the body. This makes it possible to gain much 
fuller information on the nature of the stress distribution. It also makes 
it possible to evaluate surface tractions at points other than where the 
chosen boundary conditions were applied and thus to check the quality of a 
solution obtained from the SOLVE mode and MATSOL. 
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INPUT FOR FRAC3D AND ITS AUXILIARY PROGRAMS 


Program LATTICE 

Because the input for FRAC3D is relatively extensive and suffi- 
ciently complex so that human error is easily induced, a great deal of the 
input data generation has been automated. 

The principal function of LATTICE is to generate a set of indices 
for ordering lattice points on each face of the body and for associating 
those in each pyramidal base, beginning with minimal input. It begins with 
a list of lines which the user has decided will be employed in each grid. 

Each line segment is defined by three constants: the x or y coordinate, 

and its two bounding extrema. The horizontal lines are denoted as y = C^, 

± x <_ Cy and the vertical lines are denoted as x = D^, ± y <_ D^. 

The user supplies the C's and D's in local coordinates as input for each 
face involved in the analysis. From them the program is to find the lattice 
points, enumerate them, and associate those for each pyramidal base. 

Figure 3 shows a typical grid work. Putting all the pyramid points 
first, the lattice points are enumerated in normal reading order (with 
increasing x and decreasing y) . The figure shows the indices for its case 
where space allows. The nine points associated with each pyramidal base 
are also to be taken in the same order, so that for the first pyramidal 
base in the example the indices are 210, 211, 212, 219, 1, 2, 227, 13, 14. 
Such orderings are needed in the input to FRAC3D. 

The form of the input to LATTICE is described as follows. 

Card Type 1 - Mode Card (1 card) : 

Read MODE, MFLAG 
Format (A4,70X,A6) 

MODE = SOLVE for matrix generation. Program LATTICE generates 
input data only for this mode of operation for FRAC3D. 

MFLAG = blank if no matrix print-out is desired or OUTPUT if 
matrix print-out is desired. 
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FIGURE 3. ILLUSTRATIVE LATTICE SHOWING INDICES USED BY FRAC3D 
(FRONT LATTICE FOR A/2C = 0.25, C = 0.8 a) 
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Card Type 2 - Title Card (2 cards required) : 


Read TITLE1, TITLE2 
Format 2(8A10) 

Each card contains any information the user desires to use to 
identify the run. This title also appears in subsequent runs 
of MATSOL and FRAC3D in the RESULT mode of operation. 


Card Type 3 - Symmetry Card (1 card) : 


Read JXS, JYS, JZS 
Format (3A4) 

JXS = IX1S for 
= 1X2S for 
JYS = IY1S for 
= IY2S for 
JZS = IZ1S for 
= IZ2S for 


symmetry about x 
symmetry about x 
symmetry about y 
symmetry About y 
symmetry about x 
symmetry about z 


y = 0 line 
0 plane 

0, z = z t /2 line 
0 plane 

0, z = z fc /2 line 
z fc /2 plane 


Card Type 4 - Symmetry Relation Cards (variable number of cards) : 


Read ((AXS(I,J), J - 1,3), I - 1,3) 

Read ((AYS(I,J), J = 1,3), I = 1,3) 

Read ((AZS(I,J), J = 1,3), I - 1,3) 

Format (9F5.0) 

In these arrays, the first index (I) refers to the load component 
p(I = 1), s (I = 2), or t(I = 3). The second index (J) refers to 

a pair of faces by the scheme that J = 1 for Faces 1 and 4, J = 2 

for Faces 2 and 5, and J = 3 for Faces 3 and 6. The three arrays 

are used to prescribe desired relationships between the load com- 
ponents (p,s,t) at a point (x,y,z) and the load components (p 1 , s', 
t') at a symmetrically located point (x*, y* , z'). A value + 1, 
for example, implies that p is to be equated to p', while value - 1 
implies that p is to be equated to - p' . It is important to keep 
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in mind that the loads being equated have signs determined accord- 
ing to their respective local systems of coordinates. 

As an example, suppose a slab problem is being analyzed and there 
is symmetry about the x « 0 and y = 0 planes. Across the plane 
x = 0, it is required that p = p 1 , s=“S l , t = t f , while across 
the y = 0 plane, it is required that p = p 1 , s= s f , t = - t f . 

In this case then JXS = IX2S, JYS = IY2S, and JZS is blank. Thus, 
the first card of Type 4 should contain the instructions 


J > 

* r ioo" 

AXS = ] -100 

L 1 0 0. 


on the card it appears 1,0,0, - 1 ,0,0, 1,0,0. 


The second card of Type 4 should state that 


AYS « 


I 


J > 

10 0 " 
10 0 
10 0 . 


on the card it appears 1,0, 0,1, 0,0, -1,0,0 


The third card of Type 4, for AZS, should then be omitted. 


Each array appears on one input card and there is one card for each 
type of symmetry claimed. If there is no symmetry of a given type, 
then the corresponding array card must be omitted. 

The information on the cards of Type 4, like that on the preceding 
card of Type 3 and that on the following cards of Type 5, 6, 8 and 
9 is a means for organizing the input to FRAC3D. 
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Card Type 5 (1 card): 

Read NALFI, NALFJ, NBETI, NBETJ, NGAHI, NGAMJ, NAHATI, NAHATJ, NBHATI , 
NBHATJ, NGHATI, NGHATJ 
Format (1215) 

These quantities define the crack series made available for the 
calculation. The coefficients are: 
i - 0,1, .. . ,NALFI-1, j = 

3' ., i = 0,1, . . . .NBETI-l, j = 

3 

1 = 1,2,...,NGAMI-1, j = 

o' , i = 1,2,..., NAHATI-1 , j = 

J 

1 " 0,1,... ,NBHATI-1, j = 

y’ i - 1,2,..., NGHATI-1 , j = 

J > 

Card Type 6 (1 card) : 

Read IFACE(I), I * 1,7 
Format (715) 

This card defines the coordinate systems involved in the calcula- 
tion. The value of I defines the numbered faces in Figure 1 for 
I = 1,6 while I = 7 is the crack system. If IFACE(I) = I, then 
the I— system is included in the calculation. If IFACE(I) = 0, 
then the I— system is excluded from further consideration. 

Card Type 7 (1 card) : 

Read XLO, YLO, XU, YU, ZSUBT, NU 
Format (6F10.0) 

These quantities define the geometry and material properties. 

The geometry is measured in the first face or (1) system, which 
is also treated as the global system. 


0,1, .. . ,NALFJ-1 ; 
0,1, .» • ,NBETJ-1 ; 
0,1, .. . ,NGAMJ-1 ; 
0,1, .. . ,NAHATJ-1; 
0,1, .. . ,NBHATJ-1; 
0,1, .. . ,NGHATJ-1. 
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XLO - defines Face 2 (x = XLO plane) 

YLO - defines Face 3 (y = YLO plane) 

XU - defines Face 5 (x = XU plane) 

YU - defines Face 6 (y = YU plane) 

ZSUBT - defines Face 4 (z = ZSUBT plane) 

NU - Poisson’s ratio. 

Card Type 8 (1 card) : 

Read W, RSUBC 
Format (2F10.0) 

W(~a)) = inclination of J-axis to the xy-plane. 

RSUBC = directed distance from crack circle center to the (x,y,z) 
origin, measured in the plane of the crack. RSUBC carries 
the sign opposite that of the z-component of the crack 
center coordinates. 

Card Type 9 (1 card) : 

Read XEQ, YEQ, ZEQ, XROT, YROT, ZROT 
Format (8F10.0) 

These quantities are either 1. or 0. depending on whether or not 
force equilibrium conditions in x, y, or z directions or moment 
equilibrium about the x, y, or z axes are to be imposed. 

For each face involved in the problem, the following five cards 
types must be supplied. 

Card Type 10 (1 card) : 

Read EX, EY, FX, FY 
Format (4F10.0) 

These values are the local coordinates of the points of intersec- 
tion of the crack circle and the face being described. If the 
crack does not intersect the face, a blank card must be supplied. 
The points can be entered in either order. 
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Card Type 11 (1 card for each horizontal line segment) : 

Read Cl, C2, C3 
Format (3F10.0) 

The numbers Cl, C2, and C3 are the parameters describing the line 
segment y = Cl, C2 <_ x <_ C3. The cards must be ordered in increas- 
ing values of Cl. 

Card Type 12 (1 card) : 

Read 7/8/9 EOF (end of file card) 


Card Type 13 (1 card for each vertical line segment) 

Read Dl, D2, D3 
Format (3F10.0) 

The numbers Dl, D2, and D3 are the parameters describing the 
line segment x = Dl, D2 <_ y <_ D2. The cards must be ordered in 
increasing value of Dl. 

Card Type 14 (1 card) : 

Read 7/8/9 EOF (end of file card) 

(Repeat Card Types 10-14 for each face in the problem.) 

Card Type 15 (1 card): 


Read AO, Al, A2, BO, Bl, B2, CO, Cl, C2 
Format (9F5.0) 

These numbers describe the prescribed boundary conditions on the 
crack surface as follows: 



= AO + Alp cos 0 + A2p sin 0 
= BO + Blp cos 0 + B2p sin 0 
= CO + Clp cos 0 + C2p sin 0. 
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Card Type 16 (stress points on the crack, 1 card for each value of p): 

Read p, 0^, A0 (If no crack EOF 7/8/9 card.) 

Format (3F10.0) 

Stress points on the crack are generated from this information. 

Each stress point has coordinates (p,0) where 0 successively has 

the values 0_, 0 + A0, 0 + 2A0,..., 0 , where 0 is the 

0 0 0 max max 

point of intersection between the circle with the p and the top 

surface having a positive 0. 

Note that if automatic generation of crack boundary points and 
conditions is not desired - as when an irregular array of crack points is 
to be used - then Card Types 15 and 16 should be replaced by a blank card 
in order to call for completion of the output from LATTICE. 

The punched output from LATTICE provides a nearly complete set of 
input cards for the "SOLVE" mode of FRAC3D, lacking mainly some of the 
specifications of the boundary conditions that are to be imposed. It con- 
tains first about ten cards containing titles, symmetry information, series 
structure, forces to be used, and body extent. There follows a count card 
and a lattice deck for each face to be used, a card showing a count and speci- 
fications for the crack, crack boundary condition cards and an equilibrium 
specification card. Before use in FRAC3D, each count card needs checking 
of its first field, which should show the number of boundary condition 
points for its respective face. In addition, LATTICE supplies a card with 
the coordinates of one pyramid point from each symmetry set for each face 
(except at the crack tip) to be used in preparing specifications for sur- 
face boundary conditions (after adding one near the crack tip). It also 
supplies a card showing coordinates of the midpoint of each surface rectangle 
for possible use in a checking operation after load constants have been 
found. 

The next section details the input data needed for FRAC3D and 
notes how to assign boundary-condition points and to insert counts relating 
to them. 
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Program FRAC3D 

The program FRAC3D has two modes of operation, both of which are 
used in analyzing stresses and displacements in a three-dimensional body. 
One mode of operation ("SOLVE") is used to set up a system of linear equa- 
tions which are stored or written on tape to be solved by program MATSOL. 
The second mode of operation ("RESULT") is used to generate stresses and 
displacements at selected points of the body. Input data will be described 
for both modes of operation. 

Matrix Generation ("SOLVE" Mode) 


Card Type 1 (1 card) : 

Read MODE, MFLAG 
Format (A4,70X,A6) 

MODE = SOLV - for generating the matrix 

- RESULT - for generating stresses and displacements 
(See later) 

MFLAG = blank - if no matrix print-out is desire 

** OUTPUT - if matrix is to be printed out from "SOLVE" mode. 

Card Type 2 (2 cards, always): 

Read TITLE1 , TITLE2 
Format (8A10) 

Contain descriptions of problem. User supplied for later reference. 
This information can be carried over from LATTICE. 

Card Type 3 (1 card) : 

Read JXS, JYS, JZS 
Format (3A4) 

JXS = IX1S for symmetry about x = y = 0 line 
= IX2S for symmetry about plane x = 0 
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JYS = IY1S for symmetry about y » 0, z « z t /2 H- ne 
= IY2S for symmetry about plane y = 0 
JZS = IZ1S for symmetry about x = 0, z = z t /2 line 
= IZ2S for symmetry about plane z = z^/2. 

Card Type 4 (1 card for each symmetry claimed) : 

Read ((AXS(I,J), J = 1,3), I = 1,3) 

Read ((AYS(I,J), J = 1,3), I = 1,3) 

Read ((AZS(I, J) , J = 1,3), I = 1,3) 

Format (9F5.0) for each card 

This information is normally carried over from LATTICE. A des- 
cription of its meaning is given for input to LATTICE of Card 
Type 4. 

Card Type 5 (1 card) : 

Read NALFI, NALFJ, NBETI, NBETJ, NGAMI, NGAMJ, NAHATI, NAHATJ, NBHATI, 
NBHATJ, NGHATI , NGHATJ 
Format (1215) 

These quantities are those described for Card Type 5 of the input 
to LATTICE. 

Card Type 6 (1 card) : 

Read IFACE(I), 1-1,7 
Format (715) 

These quantities are those described for Card Type 6 of the input 
to LATTICE. 

Card Type 7 (1 card) : 

Read XLO, YLO, XU, YU, ZSUBT, NU 
Format (6F10.0) 
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These quantities are those described for Card Type 6 of the input 
to LATTICE. 

-> For each value of IFACE(I), I = 1,6 (not 7) which is I (not 0), 

there is a set of cards to be read in. These are card Types 8 

and 9. 

Card Type 8 (1 card for each set of cards of Types 8 and 9) : 

Read MAXS(I), NL(I), MAXL(I) 

Format (315) 

HAXS(I) = number of stress points to be read for Face I. 

NL(I) = number of lattice points on the I— face. 

MAXL(I) = number of pyramid points on Face I. 

Note that cards of this type will be provided as output from 
LATTICE, but the counts MAXS(I) from LATTICE will need correc- 
tion to account for the number of boundary condition points 
assigned afresh by the user for FRAC3D. 

Card Type 9 (1 card for each lattice point, usually punched by LATTICE) : 

Read XL(J) , YL(J), (IN(L,J), L = 1,9) J - 1, NL(I) 

Format (2F10.0, 913) 

ttl til 

XL(J), YL(J) - coordinates of J — lattice point on the I — face 
IN(L,J) - indices of the lattice points describing the base of 
the pyramid. 

-*■ Repeat Card Types 8 and 9 for each face involved in the problem. 

Card Type 10 (1 card) (To be omitted if no crack is included.): 

Read MAXSI(7), W, RSUBC 
Format (15, 5X, 2F10.0) 

MAXSI(7) = number of stress points on the crack. 

W(~w) = inclination of y-axis to Face 1. 
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RSUBC = directed distance from crack circle center to Face 1 

origin, measured in the plane of the crack. RSUBC > 0 
if center is above first face, negative otherwise. 

There is a set of cards of Types 11 and 12 for each face in the 
problem and for the crack. 

Card Type 11 (1 card for each face in the problem) : 

Read TITLEl 
Format (8A10) 

TITLEl = stress points for Face I 

This card is provided as a spacer between the several sets of 
stress points. 

If the user augments any set of stress points, the additional 
cards may be placed anywhere within the appropriate set of cards. 
That is, the stress points are not ordered in any particular 
fashion. 


Card Type 12: 


Read XSA(J), YSA(J), SIGMA(1,J), SIGMA(2,J), SIGMA(3,J) 
Format (5F10.0) 

XSA(J), YSA(J) - stress point coordinates on I 
crack) 


th 


face (not the 


SIGMA(1,J) - value of at the stress point 

Z fi} 

SIGMA(2, J) - value of t j at the stress point 

/ • \ 

SIGMA(3,J) - value of t J^at the stress point 

In the case of the crack the coordinates to be shown are p and 9 
th 

for the I — point, and the stress components which can be pre- 
scribed are a- , T n , t , respectively. 
y ey 

Repeat Card Types 11 and 12 for each system in the problem. 
Note that boundary conditions on the body faces are usually 
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applied only at pyramid points (except for one near the crack 
tip). The coordinates XSA(J) and YSA(J) for them are provided 
among the punched output from LATTICE. 

Card Type 13 (1 card) : 

Read XEQ, YEQ, ZEQ, XROT, YROT, ZROT 
Format (6F10.0) 

XEQ = 1. if equation describing force equilibrium in the x direction 
is to be included in the system of equations. Otherwise, 
use 0. 

YEQ = same, but for equilibrium in the y direction. 

ZEQ = same, but for equilibrium in the z direction. 

XROT * 1. if the sum of the moments about the x axis are to be 
equated to zero. Otherwise, use 0. 

YROT = same, but for moments about the y axis. 

ZROT = same, but for moments about the z axis. 

The SOLVE mode of operation results in several actions. The sys- 
tem of equations and certain other needed information is written 
on TAPE 9 for use in MATSOL, which solves the system of equations. 
In addition, considerable information is written on TAPE 8 for 
use in the RESULT mode of operation. Thus, the user must identify 
two storage units which can be referenced later. 

Stress and Displacement Generation ("RESULT" Mode) 

The input data required by FRAC3D to generate stresses and dis- 
placements is listed below. The program receives a great deal of its input 
data from TAPE 8 and from TAPE 7 generated by the program MATSOL discussed 
below. 
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Card Type 1 (A4) : 

Read MODE, MFLAG 
Format (A4,70X,A6) 

MODE = RESULT, MFLAG = blank. 

Card Type 2 (215): 

Read MAXSX(l), MAXSI(7) 

MAXSI(l) - number of points in the global system at which stresses 
and displacements are to be generated. 

MAXSI(7) - number of points expressed in the crack system at 

which stresses and displacements are to be generated. 

Card Type 3 (Title card for these points, wording chosen by user) . 

Card Type 4 (3F10.0) (1 card for each point): 

Read XSA(l), YSA(I), ZSA(I) 

These are the points, expressed in the global coordinate system, 
at which results are desired. 

Card Type 5 (Title card for these cards). 

Card Type 6 (3F10.0) (1 card for each point): 

Read RHO(I) , THETA(I) , ZETA(I) 

These are the points, expressed in the cylindrical coordinate 
system, at which results are desired. 

Program MATSOL 

Program MATSOL is used after a matrix has been generated by FRAC3D 
and, together with other needed data, has been stored on TAPE 9. In its 
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ordinary version, MATSOL receives all of its input from TAPE 9 and uses it 
in solving for the load constants for all the plane surfaces used and for 
the crack if it is used. MATSOL prints out the equation residuals, the 
solution vector, and the distribution of the stress-intensity factors K^, 

Kjj, and around the crack, all normalized by division by for a 

unit load. In addition, the solution vector is written on TAPE 7 for use in 
FRAC3D operated in the RESULT mode. 

In its ordinary version, the only special consideration in perform- 
ing MATSOL is that the dimensions provided in cards 3 and 6 need to he large 
enough for the arrays A, R, C, D, B, S. All of these but A must have 
dimension as great as NT, which is the sum of the number of load constants 
which are allowed (on all the planes around the crack) plus the number of 
equilibrium conditions. A must then have dimension NT(NT + l)/2 or more. 

A special editing version of MATSOL permits special weighting of 
blocks of boundary conditions, and allows deletion of chosen load constants. 
Operation of this version requires data cards of three or four kinds as 
follow. 

Card Type 1 (1 card) : 

Read WT1, WT2 , . . . , WT 8 
Format (8FI0.0) 

These are weights to apply to equations in 8 blocks defined by 
Card Type 2. 

Card Type 2 (1 card) : 

Read IR1, IR2, » . . , IR 8 
Format (814) 

These define 8 blocks of rows in the matrix of boundary conditions, 
having weight WT1 through row IR1, then having weight WT2 through 
row IR2 , then having weight WT 3 through row IR 3 , etc. 
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Card Type 3 (1 card) : 

Read NAUX, LZ, MSR 
Format (314) 

NAUX tell how many equilibrium conditions should be enforced. 

This NAUX usually agrees with the count assigned through FRAC3D, 
but putting NAUX=0 allows deletion of equilibrium equations. 

LZ tells how many crack constants already admitted in the 
construction of the matrix by FRAC3D are to be deleted by 
MATSOL. Use of MSR > 0 permits deletion of shear boundary 
conditions on the crack by putting MSR = row number of the 
last boundary condition on a body surface (not the crack) . 

Card Type 4 (variable number of cards, none if LZ = 0): 

Read LE(I) I - 1, LZ 

Format (2014) 

This provides the column-count identification of each load 
constant to be omitted from the calculations by MATSOL. 

(This is most often applied in editing the crack function 
series. It requires advance knowledge of order numbers 
for crack constants.) 

See NASA CR 159400 for examples of why and how to edit crack 
series and sometimes surface load constants). 

In the editing version of MATSOL, card 3 should provide dimensions for A, 

R, C, D, and B as does card 3 for the ordinary MATSOL. In addition the 
sixth card of the editing version of MATSOL should provide dimension for S 
to be at least large as NT + LZ, and the seventh card should provide dimen- 
sion for JE at least as large as LZ and for JS at least as large at NT. 
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ILLUSTRATIVE CALCULATIONS FOR A SLAB WITH A SURFACE CRACK 


In order to illustrate how a stress analysis is performed with 
FRAC3D, and to show typical numbers that arise, consider a slab with a part- 
circular surface crack. For the illustration, the crack is taken to be 
perpendicular to the faces of the slab, passing into the slab to a depth 
0.4a where a is the crack radius, and penetrating 0.7 of the thickness of 
the slab. Poisson’s ratio v is taken to be 0.30. Designating the depth of 
the crack by A, the surface length of the crack by 2C, and the thickness of 
the slab by T, the defining parameters thus are A/2C = 0.25 and A/T =0.7. 

In particular, also T = 0.57142857a, and C = 0.8a. The crack radius a is 
taken as the reference length, that is a = 1. The loading is taken to be 
a uniform normal load applied to the crack surface, that is the crack- 
induced part of the loading that would arise when the slab is subjected to 
remote tension perpendicular to the crack. 

It is possible to select lattices for the front (cracked) and back 

surfaces intuitively, but it helps to consider in advance what stresses would 

arise on similarly situated planes in an infinite body having a circular 

crack under uniform normal load. The stresses o , t , and t that would 

z yz gx 

arise on those planes under a unit normal load are simply /2/tt times the 
influence functions for ^ [2,3] in rows of a matrix computable by FRAC3D 
in its SOLVE OUTPUT mode. This suggests an occasional preliminary, optional 
calculation using only minimal surface lattices (one pyramidal base on each 
surface), neglecting symmetries, and using series which need to be long 
enough only to include functions for Points where the influence func- 

tions are to be found are treated as boundary-condition points, though no 
solution for load constants is intended. From computations of this sort, 
curves of fixed levels of the influence can be interpolated, as shown in 
Figures 4 and 5*. The stresses o and t are considered there, since they 

2* ZX 

are probably the more influential ones in surface crack stress analysis. 

Front and back lattices were chosen as also shown in Figures 4 and 5, which 
hopefully should accommodate the needed freeing stresses on the surfaces. 


* Figure 4 shows one quadrant of the full lattice shown in Figure 3, the re 
complete with indices. The effective load for Figures 4 and 5 is a = a/tt/2 . 
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though further interactions between surface and crack stress systems are 
yet to be considered. 

In constructing the system of boundary condition equations, the 
surface loads were to use the lattices of Figures 4 and 5. The line seg- 
ments for the lattices are shown in the tabulation on pp 37-38. Note 
that the segments are specified for all four quadrants, not just one 
quadrant. The crack series were desired with m even, 0 <^m <^16, and 
0 < k £ 16 m, and these could be provided by putting NALFI = 17 and 
NALFJ = 9, with m-eveness and trangulation on k to be arranged later. 

Thus the calculation began by running the program LATTICE with input as 
follows. (An abbreviated coded form of this input is shown on p 39.) 


Card 

1 

2 and 3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 - 39 
40 

41 - 63 

64 

65 

66 - 82 
83 

84 - 104 

105 

106 


Data 

SOLV (in cc 1-4) 

(Titles which are shown here at the top of Page 51.) 

IX 2S (cc 1-4) and IY2S (cc 5-8) 

bbbl . bbbb . bbbb . bb-1 . bbbb . bbbb . bbbl . bbbb . bbbb . 

bbbl . bbbb . bbbb . bbbl . bbbb . bbbb . bb-1 . bbbb . bbbb . 

bbbl7bbbb9 

1 (cc5) and 4 (cc20) and 7 (cc35) 

In 6F10.0: 3.6, -3.6, 3.6, 3.6, 0.57142857, 0.30 

In 2F10.0: 0.0, 0.6 

In 8F10.0: 0.0, 0.0, 1.0, 0.0, 0.0, 0.0 

In 4F10.0: 0.0, -0.8, 0.0, 0.8 

(27 cards showing Cl, C2, C3 in 3F10.0 for front face) 
7/8/9 EOF 

(23 cards showing Dl, D2, D3 in 3F10.0 for front face) 
7/8/9 EOF 

blank card (because no crack tips are on back face) 

(17 cards showing Cl, C2, C3 in 3F10.0 for back face) 
7/8/9 EOF 

(21 cards showing D1,D2, D3 in 3F10.0 for back face) 
7/8/9 EOF 

blank card (since crack boundary conditions were to be 
special) . 
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FACE -W.-“l — — " ' 

COORDINATES OF CRACK TIPS * •Pi 1 .®! 0 . 00 ’ -0.8000000 AND. .000 0000 0.80 0 000 0 


C LINE SEGMENTS (INPUT) 



C .\ 

Cl 

C 3 

1 

- 3.6000000 

- .6 OPO 0 0 0 

3.6000000 

2 

- 2.4000000 

- 3.600003 

3.6000000 

3 

- 1.6000000 

- 1.6000000 

1.6000000 

4 

- 1*2000000 

- 3*6000000 

3.6300000 

5 

- 1 * 00000 00 

-0 .8 CO 000 0 

0.8 0000 0 0 

6 

- 0.9000000 

- 0.4000000 

0.4000000 

" 7 

- 0.8500030 

- 0 . 10 C 0000 

0.1000000 

8 

- 0. 80 000 00 

-1 *600000 0 

1.6000000 

9 

- 0*7500000 

- 0.1000000 

0.1000000 

10 

- 0.7000000 

-n .4000000 

0.4000000 

11 

- 0.6000000 

- 0 .8000 00 0 

0.8000000 

12 

- 0*4000000 

- 2.4 000 00 0 

2.40 000 0 0 

13 

- 0*2003000 

-0 .4000000 

0.4000000 

14 

*0000000 

- 3.6 000 00 0 

3.6300000 

15 

0.2000000 

- 0.4000000 

0.4000000 

16 

0.4003000 

- 2.4000000 

2.4000000 

17 

0*6000000 

- 0.3 (00000 

0.8000000 

IS 

0*7000000 

- 0.4000000 

0.4000000 

19 

0.7500000 

- 0.1000000 

0.1300000 

l.. 20 

o.aoooooo 

-1 .600000 0 

1.6300000 

21 

0.3500000 

- 0.1000000 

0.1300000 

L _ 22 

0*9000000 

- 0 .4000000 

0.4900000 

23 

1.0000000 

-0 .8 000000 

0.6303000 

l 24 

1*2000000 

- 3.6000000 

3.6309000 

25 

1.6000000 

- 1.6000000 

1.6900000 

26 

2.4000000 

- 3.6 000000 

3.6900000 

27 

" ' 3.6000000 

-3 .600000 0 

3.6000000 

OJLI NE 

SEGMENTS ( INPUT ) 




t ?' 

PZ 

93 

1 

- 3.6000000 

- 3 * 6 COO 00 0 

3.6300000 

2 “ 

- 2.4000000 

'-"3 ".6 000 00 0 

3.6000000 

3 

- 1.6000000 

- 3.6 000 00 0 

3.6900000 

' ~ 4 “ 

- 1.2000000 

- 0 *4 0 P 0 00 0 

0.4000000 

5 

- 0.8000000 

- 3*6300000 

3.6000000 

6 

“ - 0.6000000 

~- 0. 401000 0 

0. 4000000 

7 

- 0.4000000 

- 2.400000 0 

2.4000000 

8 

- 0.30 000 00 

~ : 0 . 600^00 0 “ 

0. 6000000 

9 

- 0.2000000 

- 1.6000000 

1.6000000 

10 

'- 0.1000000 

- 1.2000000 

1.2000000 

11 

- 0.0500000 

-i.oonoooo 

1.0900000 

12 " 

.0000000 

-3 .6 CO 000 0 

3.6000000 

13 

0.0500000 

-1 .0 00000 0 

1. 00 00000 

14 ‘ 

o.ioooOoo 

-‘ 1.2 COO 00 0 

1.2000000 

15 

0.2000000 

- 1.6 000000 

1.6109000 

“ 16 

' 0.3000000 

- 0.6000000 ‘ 

3.6900000 

17 

0.4000000 

-2 .4000000 

2.4000000 

18 

0.6000000 

'-0 .4000000 

0.4000000 

19 

0.8000000 

- 3.6000000 

3.6300000 

20 

1.2000000 

000 00 0 

0.4300000 

21 

1.6000000 

- 3.600000 0 

3.6000000 

22 

2.4000000 

- 3.6000000 

3.6300000 

23 

3.6000000 

- 3.6000000 

3.6000000 > 
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face ' no ; — t * 

COORDINATES O F CRACK TIPS » * 9 . 00 0 0030 , * 9.0000000 ANO * 9. 00000 0 0 , » 9. 0000 000 


C LINE SEGMENTS (INPUT) 


i 

- 3*6000000 

- 3*6000000 

3*6000000 




z 

- 2.4000000 

-■* .6 onoooo 

3.6300000 




3 

- 1.6000000 

- 1.6000000 

1.6000000 




4 

- 1.2000000 

- 1.600000 0 

3.6000000 




5 

- 0*8000000 

-1 .6000000 

1.6000000 




6 

- 0*6000000 

- 0.2000000 

0.2000000 




7 

- 0*4000000 

- 0*3000000 

0.6000000 




6 

- 0.2000000 

-0 *2 0 C n 00 0 

0.2300000 




9 

•0000000 

-“*•6 000000 

3.6100000 




10 

0*2000000 

- 0. 2000 00 0 

0.2000000 




11 

0*4000000 

• 0.8 00000 0 

0 . SO 000 00 




12 

0*6000000 

- n . 2000000 

0.2000000 




13 

0 * 8 000 00 0 

-1 .600000 0 

1.6100000 




14 

1*2000000 

- 1.6 00000 0 

3.6300000 




15 

1*6000000 

-1 .6000000 

1. 6000000 




16 

2.4000000 

-^.600000 0 

3.6000000 




17 

3.6000000 

- 3.6000000 

3*6000000 





D LINE SEGMENTS CINfUJT) 


1 

- 3*6000000 

- 3.600000 0 

3*6000000 

2 

• 2*4000000 

- 3.600000 0 

3.6000000 

3 

- 1.6000000 

-**. 6 000000 

3*6100000 

4 

- 0*8000000 

- 3.6 000000 

3*6000000 

5 

- 0 . (*000000 

- 1*6000000 

1.6000000 

6 

- 0.2000000 

-t *200000 0 

1.2000000 

7 

- 0,1500000 

- n , 2ocoooo 

0*2000000 

6 

- 0.1000000 

- 0.3 00000 0 

0*8000000 

9 

- 0. 0500000 

-0 .600000 0 

0*6000000 

10 

- 0.0250000 

- 0 *4 000 00 0 

0*4300000 

11 

•0000000 

- T . 6 00000 0 

3.6000000 

~i 2 

0*0250000 

- oVloooooo 

0*4000000 

13 

0*0500000 

-0 .6000000 

0. 6000000 

14 

0.1000000 

- 0.6000000 

0*8000000 

15 

0*1500000 

- 0.2000 000 

0*2000000 

16 

0*2000000 

- 1.2000000 

1*2000000 

17 

0 * 40000 00 

- 1.6000000 

1*6003000 

16 

0 * 8000 0 00 ** 

- 3.6000000 

3.6000000 

19 

1.6000000 

- 3 . 6 COOOOO 

3.6000000 

20 

2*4000000 

- 3.6000000 

3*6000000 

21 

3*6000000 

- 3 . 60 C 000 0 

3*6000000 


39 













40 


The major output from LATTICE is the indices for all the lattice 
points, and sets of nine indices for all the bases of pyramids. Lists of 
these indices are printed by LATTICE and again later by FRAC3D as a record 
of input it received. The list of these indices reported by LATTICE is 
shown on pp 41 to 47. 

The input for FRAC3D came mainly from LATTICE, but that was sup- 
plemented by adding five boundary condition points on the front near the 
crack tip (to be used at option) and also 80 boundary condition points on 
the crack as shown in Figure 6. Boundary tractions were assigned as needed, 
these being zero except for the normal load on the crack, which was taken 
to be unity, so that in effect the normal load on the crack became the 
reference stress. Thus the following data cards were assembled for running 
FRAC3D. (An abbreviated coded form of this input is shown on p 48.) 


Card 


Data 


1-9 


10 

11 - 315 
316 

317 - 497 

498 

499 

500 - 567 
568 

569 - 604 
605 

606 - 685 
686 


Like cards 1 - 9 of input to LATTICE (Note that in 
card 9 repunching of ZSUBT restores accuracy dropped 
by LATTICE.) 

bbb68bb305bb209 (The 68 was corrected.) 

(305 cards from LATTICE for the front face) 
bbb36bbl81bbl09 

(181 cards from LATTICE for the back face) 
bbb80, 0.0 (in cc 11-13), 0.6 (in cc 21-23) 

(Title card: STRESS POINTS FOR FACE 1) 

(68 points for boundary conditions on front) 

(Title card: STRESS POINTS FOR FACE 4) 

(36 points for boundary conditions on back) 

(Title card: STRESS POINTS FOR FACE 7) 

(80 cards for boundary conditions on crack) 

In 6F10.0: 0.0, 0.0, 1.0, 0.0, 0.0, 0.0 


It may be observed that symmetry around both the planes x = 0 and 

y = 0 was specified by Cards 4, 5 and 6. The series specified by Card 7 

called for 17 x 9 series with constants a' , , but this was intended to be 

m,k 

reduced by making FRAC3D remove the series terms for which m would be odd 
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~ LATTICE POINTS Otf FRONT FACE (PI) 


f ACi~NO. ft OUTPUT " >AXSl» 

Nil: 

30 3 HA XL 1*209 





- - 



- 






— 

- - 

-a.^oooooa 

2.4000000 

710 

211 

313 

219 

1 

i 

227 

13 

14 

-1.6000 000 

2.-000000 

211 

212 

213 

1 

2 

3 

13 

14 

15 

-0.5000000 

3.4000000 

212 

211 

214 

2 

3 

4 

233 

6 

10 

.0000000 

2.4000000 

213 

21* 

315 

3 

4 

5 

8 

10 

12 

0.5000000 

2. 4000000 

214 

315 

216 

4 

5 

6 

io 

12 

326 

l.bOOOOOO 

3. 4000000 

215 

216 

217 

5 

6 

7 

21 

22 

23 

3.4000000 

2.4000000 

216 

217 

216 

6 

7 

222 

32 

33 

230 

-o.aoooooo 

1.4000000 

2 

3 

230 

323 

3 

9 

14 

15 

16 

-o.^oooooo 

1.6000000 

3 

220 

4 

8 

9 

10 

15 

16 

18 

.0000000 

1.6000000 

23 0 

4 

231 

9 

10 

11 

16 

18 

20 

0.4000000 

1.600000Q 

1* 

2a 

5 

ir 

It 

13 

18 

30 

21 

0.5000000 

1.6000000 

231 

5 

6 

ti 

12 

236 

30 

31 

22 

-3.1*000000 

1.3000000 

219 

t 

3 

227 

13 

14 

257 

95 

56 

-1.6000000 

1.2000000 

1 

3 

3 

13 

14 

15 

249 

71 

72 

-n. aoooooo 

1.2000000 

223 

8 

9 

14 

15 

16 

239 

41 

42 

-0.4000 000 

1. 2000 000 

8 

9 

72 4 

15 

16 

17 

231 

34 

25 

-0.3000000 

1.3000000 

9 

324 

10 

16 

17 

16 

24 

35 

27 

• ooooooo 

1.3000000 

22 4 

10 

325 

17 

18 

19 

25 

27 

29 

0.3000000 

1.3000000 

10 

225 

ll 

18 

19 

20 

27 

39 

30 

0.4000 000 

1.3000000 

225 

11 

13 

19 

20 

31 

29 

30 

234 

0. 5000 000 

1.3000000 

11 

13 

326 

30 

21 

23 

50 

51 

240 

l.bOOOOOO 

1.3000000 

5 

6 

7 

21 

22 

23 

84 

85 

254 

3.40 00000 

1.2000000 

6 

7" 

222 

33 

23 

330 

114 

115 

256 

-0.4000000 

1.0000000 

15 

16 

17 

231 

24 

25 

41 

42 

43 

-0.3000000 

l. ooooooo 

16 

17 

228 

34 

25 

36 

235 

31 

32 

-0.1000000 

1.0000000 

17 

239 

18 

25 

26 

37 

31 

32 

34 

•ooooooo 

1 . ooooooo 

22 6 

13 

229 

26 

37 

28 

32 

34 

36 

0.1 000000 

1.0000000 

18 

339 

19 

27 

28 

29 

34 

36 

37 

0.3000000 

1.0000000 

22 9 

19 

20 

28 

29 

30 

36 

37 

236 

0.4000000 

l.OOOOOOC 

19 

20 

21 

29 

30 

334 

49 

50 

51 

-0.3000000 

3.9000030 

24 

25 

26 

2 35 

31 

32 

42 

43 

44 

-O.lOQOOQQ 

0.9003000 

25 

36 

233 

31 

32 

33 

43 

44 

45 

-0.0500000 

0.9000000 

26- 

233 

27 

32 

*3 

34 

237 

38 

39 

•ooooooo 

0.9000000 

212 

27 

23 3 

33 

34 

35 

36 

39 

40 

0. 0500000 

0.900000ff 

27 

233 

38 

34 

35 

36 

39 

40 

238 

0.1000000 

0.9000000 

733 

23 

39 

35 

36 

37 

47 

48 

49 

0.2000000 

0.9000000 

26 

29 

30 

36 

37 

236 

48 

49 

50 

-0.0500000 

0. 8500000 

33 

33 

34 

237 

38* 

39 

44 

45 

46 

•ooooooo 

" 0.8500000'"* 

33 

34 

35 

38 

39 

40 

" 45 

46 

47 

0.0500000 

0.8500000 

34 

35 

36 

39 

40 

238 

46 

47 

46 

-0.50 00 030 

* 0.8000000* 

14 

15 

16 

239 

41 

42 

71 

72 

73 

-0.4000000 

0.8000000 

231 

24 

25 

41 

42 

.43 

245 

62 

63 

-0.3000030 

0.8000000 

23 5 

31 

32 

42 

43 

44 

243 

55 

56 

-0.1000000 

0.8000000 

31 

32 

33 

43 

44 

45 

55 

56 

57 

-0.0500000 

0.8000000 

237" 

"38 

39 

44 

45* * 

" 46 

“ 241 

52 

53 
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- 0. 2000000 
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148 
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0 . 1500 000 

- 0.2000000 
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149 
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0.2000 000 

- 0 • 2000 000 
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150 
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- 0,8000000 " 
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~~ 0 " 
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- 0.3250 000 

- 0.4000000 
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152 
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LATTICE POINTS ON BACK FACE (P3) 

0.0250 000 - 0 * 4000000 0 0 fi 0 153 "" " 0 0 0 

0*5000 000 -0.4000000 0000 154 0 0 0 

-0.20 00 000 -0.6 00 0000 O'” "0 0 ” 0 " 155 0 0 0 

-0.05 00 000 -0.6000000 0 0 0 0 156 0 0 0 

0.0500 000 -0.6000000 0 0 0 “ 0 ' 157 0 0 0 

0.2000 000 -0. 6000000 0 0 0 0 158 0 0 0 

-1.60 00 000” -a", ft 00 0 0 00 o o 0 0 159 “O' o “o' 

-0.1000 000 -0.0000000 0 0 0 0 160 0 0 0 

o.ioooooo -a.aoooooo o' o' o'" o i6i o o o 

1.6000000 -O.ftOOOOOO 0 0 0 0 162 0 0 0 

-3.6000 000 -1.2000000 ' 0 0 Q 0 163 '0 0 ' 0 

-0.2000000 -1.2000000 0000 164 000 

’ 0*20 00 000 “”-1.2000000 0 3 ff 0 165 0 ~<f 0* 

3.6000000 -1.2000000 0000 166 000 

-1.6000000 -1.6 00 0 000 ' 0 * '0 0 0 T 167 '0 0 '~0‘ 

-0.40 00000 -1.6000000 0 _ 0 0 0 168 0 0 0 

, 0. 40 UO 000 -1.6000000 0 0 “ 0 0 169 0 0 0 

’ 1.6000 000 -1.6000000 _ 0 0 __ 0 _ „ 0 __ 170 0 0 _ 0 

-3.5000000 -2.4000000 0 0 0 0 171 0 0 0 

3.6000 000 -2.40 0 0 000 0 0 0 0 172 0 0 0 

-3.6000000' -3.6000000 ~0 0 0 0 173 0 0 0 

-2.4000000 -3.6000000 0 0 _ 0 „ . __ 0 174 0 0 0 

-1.6000000 -3.6000000 0 0 0 0 175 0 0 0 

-0.8000000 -3.6000 000 0 _ 0 0 0 176 _ 0 _ 3 0 

.0 000000 -3.6000000 0 0 0 0 177 0 0 0 

0.8000000 -3.6000000 0000 178 0 00 

1.60 00 000 -3.6000000 0 0 Q V * 179"' 0 0 ' 6 

2.4000 000 -3.6000000 _ 0 0 _0 __ 0 _ 180 0 0 0 

3.6000 000 -3.6000000 ~ 0 * 'O 0* 0 181~ 0 0 0 




FIGURE 6. BOUNDARY CONDITION POINTS AND CHECK POINTS 
FOR SURFACE CRACK WITH A/2C =0.25 
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or for which m + 2k would be greater than twice the maximum k (that is 
greater than 16). To accomplish this removal, 18 special cards (as have 
been provided with FRAC3D) were inserted into the subroutine AXYZ, between 
the instructions "100 CONTINUE" and "150 CONTINUE". The introduction of the 
five boundary points at or near the crack tip was made with the understand- 
ing that this would offer a choice among possible boundary conditions to be 
imposed near the tip, the boundary points that were used can be found hand- 
written as a table given later (Pages 53 to 56 ) . Card 686 called for 
satisfaction of equilibrium in the z direction among the freeing loads 
applied to the front and back faces. 

The allowance made for symmetry and the deletions of the unwanted 
crack constants are shown on Pages -50 through .52 as printed by FRAC3D, 
using the system of notation described earlier. The beginning of this chain 
shows, for example, that influence functions for the first load constant 
were added to those of the nineteenth, and those for the second were sub- 
tracted from those for the twentieth. The zero entries show which load 
constants were retained for solution and counting shows that 303 of the 
possible 1404 were retained. The other standard printed output from 
FRAC3D consists mainly of a listing of the boundary condition points that 
were used showing their coordinates used in calling the crack and surface 
subroutines in the program. 

The major output from FRAC3D is the matrix to be solved, but since 
it is typically large, it is simply placed in a computer storage unit iden- 
tified as TAPE 9, including the symmetry chain. The input which led to the 
matrix is stored also in a unit called TAPE 8. The solution for the 
implied system of equations is obtained then by use of the program MATSOL, 
which calls for the information in TAPE 9. After the system is solved, 
MATSOL stores the load constants on a unit called TAPE 7 and it prints 
output of several types. First, it prints the residual from each of the 
equations used in the least square process. For the illustrative example, 
these residuals are shown on Pages 53 through 56, showing fits for a^, 

t , x , or o', x' , t' or , t , t at each point indicated for fit- 
yz* zx* z’ yz zx J fir p r 

ting. Note that these residuals are affected by the weights that were 
applied because of symmetries specified. It may be seen that altogether 



SYMMETRY CHAIN (PI) 
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562 equations were used, but where the residuals show E-ll, E-12, or E-13, 
the corresponding equation involved only essentially vanishing influence 
functions. The many residuals showing E-Q2 and E-03, being << 1, indicate 
good fitting at these points. These residuals provide a means for studying 
the quality of a solution, but they do not tell all that might be desired 
since satisfaction of boundary conditions is desired everywhere on the front, 
back and crack surfaces. 

In this example, the editing version of MATSOL was used to select 
boundary conditions to be fitted near the crack tip. The weights of the 
various rows of the matrix (not of the printout) are shown near the begin- 
ning of the printout, and these show that here rows 1 to 4 and 8 to 15 were 
assigned weight zero, so that rows 5 to 7 provided the boundary conditions 
used near the crack tip. The occasional large residuals among the non- 
fitted rows there suggest the difficulty of getting good fitting all around 
the crack tip. 

MATSOL uses the auxiliary information provided by FRAC3D to 

identify the load constants, and the identified forms are printed as shown 

on pages 58, 59, 60. Then it uses the crack constants to compute stress- 

intensity factors for the three standard modes at many points along the 

crack front as shown on Page 61. In the illustrative example, 40 intervals 

were used from the crack root to the crack because FRAC3D contained a 

statement NO = 40 which was carried forward to MATSOL. The highest value 

of 0, of course, is arccos (r /a), which here is arccos 0.6. 

c 

In this illustrative calculation, the computer proceeded directly 
from performance of FRAC3D to performance of MATSOL, Then it returned 
directly to FRAC3D with data calling for operation in the RESULT mode in 
order to evaluate stresses at a set of points specified in the global 
coordinates and a set specified in the crack coordinates. The data pro- 
vided for the calculation in this latter mode were 

Card Data 

RESULT (in cc 1-6) 

bbl00bbb71 (counts of points in global and crack systems) 
(A title card: STRESSES AND DISPLACEMENTS FOR POINTS 
ON FACE1) 


1 

2 

3 
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PYR*H10«L LOOPS OH F«CE NO» 1 


GLOBAL 

COORDINATES 



SURFACE LOADS 


X 

y 

z 

p 

S 

T " 

0. LOGO 

2.4000 

Q.QuOG 

.2&95966622E-01 


-• 19916825 COE -02 

• 8C09 

2. LOGO 

c.ooco 

.317L59LL29E-01 

-. 34 30 4307 31 E -02 

•44360376656-02 

1.6DU0 

2.L000 

C.0300 

.33148265236-01 

•28115C5140E-02 

• 1Q669235 45E-0 1 

2 • LI 00 

2. L COO 

0. OOGO 

•22972973156-01 

• 757569 86 08E -0 2 

•B169944705E-CZ 

u • CO 03 

i.eo:c 

C.OLGG 

• 12958 12 3526- j 1 


-•3534425344E-C1 

• 4CG0 

1.6000 

0 • G w C 0 

•298L739691E-01 

-.21G1223136E-01 

-. 27025693956-0^ 

« 6 u J J 

1*6030 

C. 0 jUO 

•35730 1387 7£ -01 

-.168076088 9E -01 

.82736210616-02 

0. .COO 

1 • 20 GO 

.c.ccoo 

-*4G4 47 316826-01 


-.97354454 75 6-01 

.20 (iu 

i.2w00 

0 • 0 u G G 

• L10U29G 1836-0 1 

-.40C54366186-01 

-.9033 37 36 42 £-01 

• L w u 3 

1.2GlC 

G.ObuO 

.51186423886-01 

-.7697936 29 4E-01 

-.10793838746-01 

• 60 GO 

1*2300 

0.0000 

.18866421376-01 

-.5537870 940 E-C 2 

.56812554736-01 

1.6LL3 

1*2000 

G. 0000 

.47644612506-01 

• 300 230 6660 6-01 

•2832342241E-01 

2. LOGO 

1.2GJG 

O.OOCO 

•39362641 186-01 

.23724723076-01 

•1J 729467546-0 1 

0.0000 

1.0000 

0.0000 

-.73136767916-02 


-.23561999946+01 

• 10 Cu 

l.uOCO 

0.G3C0 

•79841 94 974E-Q1 

-.49394558056-01 

-.170G8158 536+00 

• 20 00 

l.OOGC 

C.OCOO 

•996766 868 46 -Cl 

-.1C57655225E+G0 

-.23928254166-01 

• LOCO 

1.0030 

C. 0 GOO 

.40794583496-02 

-.53813371726-01 

.11593393176+00 

o.coco 

.9000 

C.GGJG 

• 22t 3u 0 1 0 66E *0 0 


-.41115287746+00 

• C50Q 

• 9UU0 

o.oauo 

• 22130 J5 2316 + OG 

-.B2eoi5i279E-01 

-.29469353386 +0 0 

• 1003 

.9000 

c.ooco 

.16498 312476+00 

-.10809597776+00 

-.10185239036+00 

• 2b 00 

. 90C0 

0.03GC 

.18293741626-01 

-.5740 60 2 336 E-C 1 

♦12D69G89 896+00 

u« CO 00 

• 8530 

o.ocuc 
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-• 5 7532323 60 h +00 

« i 5 0 w 

.650 3 

T. GGOG 
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G . 0 bO G 
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• 8000 

C. 0 JCU 
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. Lu C3 

• 6 3 C 0 
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. 6u C 0 

• 8 G G 0 

0 • 0 uO 0 
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y&m ** §1 IliltlLltl Jll 
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• CEO] 

• 7 J JO " 
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• 7 C 00 
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0 . ulCJ 
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PTRAHIOAl LOADS ON FACE MO . 4 








GLOBAL 

COORDINATES 



SURFACE LOAOS 


X 

Y 

z 

P 

S 

T 

o.ucco 
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.025 0 
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Card 
4 - 403 

404 

405-475 


Data 

(100 cards specifying x, y,z in 3F10.0 to be used in 
finding stresses) 

(A title card: STRESSES AND DISPLACEMENTS FOR POINTS ON 
FACE 7) 

(71 cards specifying p,0,^ in 3F10.0 to be used in 
finding stresses) 


The output from the calculation in the RESULT mode is shown on 

pp 63 to 66. It may be noted that points listed as Face 1 are really the 

midpoints at all the rectangles on Faces 1 and 4, so that they constitute 

a fairly complete set of check points for the fitting on the plane faces. 

The columns providing the check are those for o , x and x , and inspec- 

z yz zx 

tion shows that all those entries are a small percent of except in some 

small rectangles close to the crack tip. The points listed on Face 7 are 

the chosen check points on the crack. Here should be near unity, while 

(misprinted as TAUTR) should be small compared to unity. 

t these conditions too are well satisfied, so the 

analysis provided here is taken to be acceptable. The values of W listed 

for Face 7 provide a map of the crack opening displacement after they are 

multiplied by 2ya/a , where y is the shear modulus, and a is the applied 

o o 

uniform crack load.* 

The illustrative case described here is discussed also in Reference 2, 
as the accepted analysis for a normally loaded crack with A/2C = 0.25 and 
A/T = 0. Much of the tabular information here is shown there in graphical 
form. 


x (TAUTZ) and x 
Inspection shows t 


* The program for getting the other displacement components still needs 
checking for consistency while combining contributions for surface 
constants with those from crack constants. 
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COMMENTS ON USE OF FRAC3D 


The report to which this manual is a supplement discusses several 
aspects of the accuracy of the program FRAC3D. It is observed there that 
the formulas underlying the program seem dependable, as do the parts of the 
program that have been used for finding load constants, stress intensity 
factors, and crack opening displacements for plates or bars with cracks 
perpendicular to the front face and loaded normally or obliquely. The 
program needs checking where it combined displacements from crack and sur- 
face loads. It may be possible that some errors persist in still unused 
parts, such as those relating to obliquity of the crack, but the program 
probably must be used in new kinds of calculations in order to reveal any 
such errors. 

The time consumed by a calculation depends greatly on the amount 
of detail put into the calculation. It is possible, however, to estimate 
the central processor time required by considering LATTICE, FRAC3D, and 
MATSOL separately. The part done by LATTICE usually is brief. Thus, using 
the CDC Cyber 73 at Battelle's Columbus Laboratories, a run of LATTICE to 
get the latticework shown in the above illustrative example would require 
about 23 CP seconds. The time used by FRAC3D varies roughly in proportion 
to the product of the total number of lattice points (for both faces) times 
the number of boundary condition points used. In the illustrative example, 
there were 486 lattice points (that is 305 + 181) , and there were 175 
boundary condition points (that is 68 for the front, 36 for the back and 71 
for the crack) , and the time used by FRAC3D in assembling the matrix was 
700 CP seconds. The time required by MATSOL varies roughly as the cube of 
the number of unknown load constants. In the illustrative example there 
were 303 load constants, and solving for them required 455 CP seconds. It 
may be added that there are probably several improvements which could speed 
these calculations substantially, but they have not been pursued yet because 
the main thrust of the investigation has been related to getting depend- 
ability in the results rather than speed. 
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Attainment of dependability in results computed by FRAC3D is 
closely tied to the way the user designs the calculations to be made. Some 
of the elements of good design are built into the illustrative example 
provided here, such as deletion of odd values of m when treating a crack 
that is less than semicircular, and putting emphasis on the latticework on 
the front face. (Even when back surface effects were being sought, the 
design of front-face latticework seemed more critical than that for the back 
face.) The many cases discussed in the accompanying report [2] illustrate 
more principles of good design for the calculations and also show the 
importance of extensive checking of non-fit ted boundary conditions, since 
unchecked results can easily be poor. It is not to be thought that this 
situation reflects any unique weakness of FRAC3D; it is quite likely that 
parallel problems exist in most calculations in which solutions of two 
basic kinds (e.g. crack and surface analyses) are merged by fitting at 
discrete points. Problems of this sort probably characterize the current 
state of the art in fracture analysis. 
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APPENDIX A 

THEORETICAL BACKGROUND OF FRAC3D 


Overall Plan for Analysis 

The purpose for Program FRAC3D is to analyse three-dimensional 
stress fields around a crack with a circular front situated in a body which 
may have up to six plane faces. The crack may be wholly or only partly within 
the body, so that surface cracks may be treated, and much flexibility is to be 
allowed in the orientation of the crack and in the kind of loads applied to it 
and to the faces of the body. 

The analysis is based on a merger of one for stresses around are 
arbitrarily loaded circular crack in an infinite body, and on one for stresses 
due to arbitrarily distributed loads on the surface of a half space. Analyses 
of these two kinds, described in terms of arbitrary load constants, are to be 
merged by choosing load constants so that appropriate boundary conditions are 
satisfied on the relevant crack and body surfaces. The merger is to be accom- 
plished by the boundary point least squares technique. Thus one goal in using 
FRAC3D is to assemble equations for imposing boundary conditions at chosen 
points sampling all the appropriate surfaces. The least squares solution of 
these equations is to be performed by the auxiliary program MATSOL, which 
solves for the load constants and uses them at once to compute stress inten- 
sity factors. The evaluation of stress components at any point closely 
resembles the construction of the boundary condition equations, so it is to 
be accomplished by returning to FRAC3D after MATSOL has evaluated the appro- 
priate load constants. 

A discussion of the principle equations used in these calculations 
follows. Fuller explanation of their derivations are given in References 
1 and 4. 


Stresses Around a Loaded Circular Crack in an Infinite Body 


Since it is planned that the stresses in a cracked, finite body will 
be analyzed by merger of separate analyses of stresses from crack loads and 
from body-surface loads, one may begin by considering a crack in an infinite 
body, A form of crack which is particularly desirable for study is the 
circular crack, which is genuinely three-dimensional and for which much 
analysis is possible because of the abundance of mathematical theory relating 
to geometries that are basically circular. Therefore, this research began by 
considering stresses around a circular (penny-shaped) crack in an infinite 
body. The loads to be applied to it in the analysis are those which the body 
would transmit there if there were no c rack, but which can not be transmitted 
across the crack, so that they are in effect the crack’s addition to the loads 
on the body. In consideration of the eventual merger with effects from fairly 
arbitrary effects from body-surface loads, it is desirable to consider crack 
loads which are as general as possible. Therefore, the crack loads to be 
considered here embrace all that can be transmitted across a circular, plane 
section of a body, namely arbitrary distributions of normal loads and tangen- 
tial loads in two directions. 

A beginning for the present analysis was found in Muki's treatment 
[5] of stresses in a half space subjected to normal and tangential forces 
arbitrarily distributed over a circular section of its surface. Proceeding 
from this foundation, two half spaces were considered, being joined continu- 
ously everywhere except over a circular section of their interface, and with 
equal distributions of three load components applied to each half space on 
that circular section. Since Muki’s analysis was performed by expressing the 
load distributions by Fourier series and manipulating conditions implied for 
each term by use of Hankel transformations, this new analysis soon required 
the solution of many systems of integral equations with Bessel functions in- 
side the integrals. Fortunately it was found possible to solve these systems 
analytically, as is shown References 1 and 4. The final formulas found for 
the stresses and displacements are reasonably compact, though they do make 
free use of integrals of a seemingly formidable type. The discussion that 
follows first shows the stress and displacement formulas, and then it shows 
means to make those integrals efficiently calculable. 


Formulas for Stresses Associated 
with the Crack 


The analysis outlined above showed that all the stresses and 
displacements around an arbitrarily loaded circular crack in an infinite 
body can be expressed as simple combinations of integrals of the form 

00 

* J e ^ J M* p ^ J N+l/2^ d? • 

0 

Here cylindrical coordinates (r,0,|) are used, the crack lies in the region 
where r < a and ^ * 0, and p * r/a, £ * j/a. In order to organize the 
expressions that have been found for the stresses and displacements around 
the circular crack for computing, it is desirable to introduce several sets 
of auxiliary functions. In doing this, allowance for the appearance of 
negative as well as positive values of the axial coordinate £ is made through 
use of the notation 
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Letting v be Poisson's ratio, the first set of auxiliary functions is 
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Then the following functions appear in the expressions for the displacements 
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The following functions appear in the expressions for stresses. 
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Using the above functions of p and C, and letting a be the crack radius and 
U the shear modulus*, the expressions for displacements and stresses given 
in the cited report, after extension to include antisymmetric as well as 
symmetric loads, become: 
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* Note that p = E/[2(l + v)], where E is Young's modulus and v is Poisson's 
ratio. 
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It may be noted that several indeterminacies may arise in the 
Equations (5) for points along the axis where p « 0, namely in the ratios 
Um k, i^P and v m,k, i^P* These indeterminacies can be resolved readily by 
1 'Hospital's rule, using the relation. 



which follows ftom a recursion formula for Bessel functions. In particular, 
one finds that 
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These expressions can replace the indeterminate expressions in the Equations 
(5) for p = 0. 

In the Equations (6), the coefficients a ? . , 3 ? . , y 1 , , a 1 , , 

/s A m j k m , k. m , k. m , k. 

3 T , , and y T - are weighted sums of earlier constants first introduced as 
m,k m,k 

coefficients of Fourier-Bessel expansions for the load distributions on the 


crack. Those earlier constants (called a 


m,n 


3 , and so forth) were 

m,n 


arbitrary since they were introduced to describe arbitrary load distributions, 

so the newer constants (a 1 , , 3 r . , and so forth) are also almost all 

m,k. m,k 

arbitrary. An exception to arbitrariness for them is a mathematical require- 

* ^ t 

ment to take Y0>k = Y0>k = 0 for all k, these being constants that arise as 
fictions undetermined by the loads on the crack. The remaining newer con- 
stants can be evaluated in various ways, and in particular it is proper to 
seek them by assigning values to stresses or displacements on chosen points 
anywhere in the body and then solving equations such as those in Equations 
(6) for the desired constants. This is essentially what will be done in the 
hybrid stress analysis except that there contributions from more than one 
loading system will be considered simultaneously. A special point to observe 
is that all the functions multiplying a’ in the Equations (6) vanish 

U ^ K 

identically for all values of k, so that these constants cannot be found 

from specified stresses or displacements. The reason for this is that these 
^ . 

constants, , are also mathematical fictions, undetermined by loads on the 

* V ^ V A 1 

crack. Apart from the yo,k 5 a 0,k> and Y0 k> all th e coefficients a ? , , 3 1 , , 

, y^ k are physically meaningful and may be used in analyses of stresses 

of bodies having circular or part-circular cracks. 
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Knowledge of the coefficients a 1 , , 3’ , , , y* , can lead at 

in K in , K. 

once to an evaluation of the stress intensity factors of all three modes* 
In particular, if the stress intensity factor of the first mode is defined 
to be 


1 im 


k i B V r > 9 > 0) > 


( 10 ) 


then (as was shown in [1] and [4]) 


kj = ,/a. ^ (<% cos m 8 - O' £ sin m 0 ) , ( 11 ) 


m=0 


where 


< - fl l - -/! I • (i2> 

k^O k-0 

The stress intensity factors of the second and third modes are defined to be 


lim 


k 2 - T ^r > 


lim 


^ - r-*a+ (r,e,°) , 


(13) 


and for them it was found that 


k 2 = [(p" - y" ) cos me - (§ " - Ym ) sin me] , 


m=0 

CD 


k 3 «= v/a £ [(3m + Vm ) sinm6 + (e"+Ym ) cos me] , 


m=0 


where 


(14) 


K " I ( ' 1 )k 0 m > k ‘ ^(2-v) ’ Y m = ^I ( " 1)k Vk * (15) 

k=0 


k«=0 
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with correponding definitions for g" and y" starting from g 1 , and y* . 

m in m , tc m , k 

Knowing the coefficients in Equations (11) and (14) allows calculation of 
stress intensity factors for as many values of 0 as the analyst desires. 
Reference 4 provides interpretation for individual coefficients of those 
kinds. 

Although the major calculations contemplated here involved simulta- 
neous use of several or many load constants, it is helpful to know the effects 
associated with some particular constants that correspond to particular 
simple crack loads. Among these are a* n (~ uniform normal load), a f n 

V j U X y U 

tilted normal load), yj ^ (~ uniform, unidirectional shear) and ^ 

(~ uniform twist). These cases are discussed in Reference 3. 

Reformulation of the Crack-Function Integrals 

Although the series in the Equations (6) may not always be long, 

n* 

there are nevertheless many of the integrals ^(p,C) required since many 
combinations of the indices n 1 , M, and N appear there. Moreover, one should 
not be unduly limited in how large m and k may be or in how many p T s and 0 f s 
may be used. Thus, the workability of this analysis depends on how readily 
the integrals ijj ^(p,£) can be evaluated. It is possible with modern 
computers to perform the needed integration directly, but the time to do this 
becomes prohibitive because of the number of integrals needed. . To overcome 
this problem, reformulations of the integrals into more readily computable 
forms will be given. 

Some help for reducing the number of integrations needed can be 
had from two recursion formulas which follow directly from familiar recursion 
properties of Bessel functions. These recursion formulas are 



2(Mfl) 

P 





(16) 


n'-l 


M 


ilJ+ 2 (p,C) = (2N+3) ^(p.C) - I M>N (P,C) 


(17) 
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Use of these formulas makes it easy to compute a large set of functionals for 
a given p and £ if only a suitable base set is computed first, such as all 
those with H = 0 or 1 and with N up to some suitable value. For N > 1 it is 
possible also to omit the integrals for n' = 1 or 2, since those integrals 
too can be found by recursion. 

Values Where £ > 0 . To get efficient means for computing the base 
set of integrals for a given p and £, the integrals were recast through use 
of integral representations for the Bessel functions J^( p§ ) and 
in Equation (1). After considerable manipulation (Reference [1]), it was 
found that if £ > 0, the integrals for n f = 0 and M = 0 or 1 can be rewritten 
as 


jO , D r ) 

I 0,N (P * C) 4^fp 


J R(x,p,0‘ 


J&- p) 2 +C 2 N 


• P.,(x)dx + 


(18) 


. r R(x . r) 

4^p£ J ( ,P,C) 

-1 


x+p 


«/(*+ P ) 


• V x)dx . 

+c 


ij N (p,C) - < -1 >~ -1 J {jl pC . R(x,p,C) • 4 + - X( ' 2 + P - ) 2 - } • P N 00dx + 

* 4 /c .i 1 y<x+ P r+r J 


(19) 


. r-ir^iw-n 1 r , 1 „ , 

4 ^» ! l * (x -°- 0 'jzzw - ■ 


where 


R(x,p,C) 



V 


) + A^v) 


2 ,, 2 2 
+ 4 p c 


( 20 ) 


P N (x) is a Legendre polynomial of the form 


,N 


V x) 


2 N Nl dx N 


[(^)1 ■ i 


r*=0 


( 21 ) 
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and 



if N is even 


if N is odd 


( 22 ) 


These formulas are well suited for computation since their sensi- 
tive features are controllable. Thus, though many significant digits may be 
lost by subtraction in computing the polynomials P^Cx), those functions are 
free of p and ? so they can be precomputed carefully without repetitious 
expense. The fluctuations of the P^(x) also waste significant digits during 
the integration, but the polynomial nature of the fluctuations provides 
compensation by making them conform well to the approximations of the inte- 
grand used in the standard Gauss integration if its order is high enough. To 
complete the base set of integrals for £ > 0, one may use the following 
integrated forms for the needed cases having n T = 1 or 2. 





(23a) 



(p,C> 



(p,c> 



(23b) 


l}, 0 (P,C) 

«io<p.O 


' js; t c[ ( 1 “ l ) ' 

^.o'P.O - + ^=^5 [(c 2 +i+p) + (c 2 «-p)=3 • 

■ - 2 — 6}=5( i+c DT a+p)+a - c)c2 J +4b i[ <1+p)4<1 - p,c L]} • 

■ • <■ -a 


(23c) 
(23d) 
(23e) 
(23f ) 
(23g) 



A-13 


2 

I? 1<P.C) • - - 7=4 — 6 {'!(> -«3[atp>+<i-p)«3 + 4»i[a + p)-«-p)' 6 J} + 

4^2np C b ! c i 


(23h) 


A Tt P 0 ! 


H(‘-3 


Here 


V) ♦ /(i-p V) 2 ♦ 4p 2 c 2 '' 
11 mJ 2 [a-p> 2 + c 2 ] 





(24a) 


(24b) 


(24c) 


The formulas (18) through (24c) provide a suitable base for computing the 
needed integrals for £ > 0 (that is off the plane of the crack), and use of 
the recursion formulas (16) and (17) can complete their computation. The 
integrals needed at places where £ < 0 are to be computed in the same way 
since £ appears inside the integral only as | £ I . 

The evaluation method for the crack function integrals as described 
above is limited not only in requiring £ > 0, it also fails if p = 0 since 
several of its expressions are then indeterminate, and in Equation (16) it 
involves an explosive recursion if p is sufficiently small. Cases with £ > 0 
and p = 0 are tractable, however. Thus one may note first that the original 
integral for N (P»?) vanishes for p = 0 if M > 0. If M = 0, then evalu- 
ating the limit R(x,p,£)/p as p -*■ 0 in Equation (18) shows 
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m x 

(-1) L2J [14-(-1) n 1c j 


2 J 2 n 


vi> 

2.-2 


dx - 


(-D 


ti-(-i) N ] r xP n ( x) 


* +c 


2 /sr 


i 


2,-2 


dx 


Ji * +e 


(25) 


and values computed from this relation provide the familiar basis for recur- 
sion to get all the required values of Iq ^(0,£;) provided one includes the 
following additional limiting values for p 0. 



_ ^ 2 /tt 

i+e 2 


(26a) 


i o 1 .i (o * £) -yi( 


Z ( arctan ** 


TeT 


-kL) 

i+ c 2 


Vo <°'0 


2 g y^/rr 

(i+c 2 ) 2 


(26b) 


(26c) 



-2^ 

(1+C 2 ) 2 


(26d) 


If C > 0 and p is not zero but small, then it seems natural to 
seek evaluations made by modifying those for p = 0, and indeed such evalua- 
tions are possible. Thus one can show that 



” (-l) s /P\»2s „ 

* (M^2s)l \2/ M+2s s 


n'+Mf2s 

0,N 


( 0,0 • 


Use of such a p-series involves evaluation of integrals with p - 0 and n f 
increasingly large. Integrals with large n f have potential difficulties 
with convergence, but those difficulties can be resolved so that series of 
this form are useful. The forthcoming paper on evalution of the integrals*, 
treats these matters in detail, and the methodology described there has been 
included as an alternate procedure in FRAC3D so that integrals with small 
positive p and C > 0 can be computed efficiently. The reader is referred to 
the forthcoming paper for details. 


* J. C. Bell, "Evaluation of Integrals Involving Products of Bessel Functions 

Having Application to Crack Stress Analysis". 
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Values Where £ = 0 or is Small . For points on the plane of the 
crack (for which £ = 0), special formulas are needed, and their forms depend 
on whether the point is on the crack (with p < 1) or beyond it (with p > 1) . 
One simplification available is that all the integrals with n* = 2 are multi- 
plied by C in the formulas (5) and (6) for the stresses and displacements, so 
that they can be shovm to need no immediate consideration. 

With this in view, a suitable base of integrals for starting the recursion 
for t; = 0 and p > 1 can be shown [1] to be 


/ 2 . 1 
f “ 7 f= arcs in — 

p 

0 


X O(N <0> 0 ' 

(p > 1) 




N 

. 2 , 


for N ■= 0 , 
for N odd , 

1 N-l N 


N+l 


(27) 


i a> . r t x x dt 

— KU.1 « ' ' 

Jhi 2 


« , , N+l 

1 ■ f 0 


for N even and > 2 


(P > 1) 



for N *= 0 , 
for N «= 1 , 


for N even and > 2 , 




(-1) 2 (1-3-5-...-N) p 

2^ /2 (^)l p N+1 J 0 


N-2 N+l _ N+2 

t 2 (1-t) 2 (l-tp' 2 ) 2 dt 


(28) 


for N odd and > 3 


j 
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using also these four special cases with n 1 = 1, 


(p,0) «= 0 , 


(p > 1) 

(29a) 

2 r i 

(p,°) = “7 t=i arcsin - - 
t L p 


(p > 1) 

(29b) 

(p,°) = nr-y , 

^/2n p 7p -1 


(p > 1) 

(29c) 

(p,0) = 0 


(p > 1) 

(29d) 


For points on the crack surface itself (with £ *= 0 and p < 1), 
most of the crack functionals become either polynomials in p or polynomials 


multiplied by . In order to express coefficients of the polynomials, it 

is helpful to use the Pochhammer notation (cy)^: 

(oO n = O'(cH-l) (o+2) • . . . • (cH-n-1) , (oOq m ^ • (30) 

For £ « 0 and p < 1, a suitable base set of functionals for starting the 
recursion can then be derived [1] to be 


x 0 yp.°> 

<p < i) 


1 0 s ?) (-f ) . 2 " 

[1*3*5 It t V l _n_ n 

^(N+1)2 N/2 (|): ^ (1) n n! 

2 (Wl)/2^, ¥ (if) (if*) P 2 ” 

JZn (1*3*5* *N) L (1) n! 

a n 


for N even , 


for N odd 
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I 1,N<* 5 ’ 0 ^ 

(p < 1) 


f jk p [ L • 


N-2 


mm .' 1 

[ 1 • 3 - 5 •(N-1)]Zj (2) nl 


for N = 0 , 


for N «= 2,4,6,--- , 


n«=0 
N-l 
2 


(32) 


(1*3*5* *N) 




Q_np_ y 

N=i'ii L> 


n*=0 


(¥) (¥) p 2 " 

n n 

< 2 V : 


for N odd , 


using also these four special cases with n 1 = 1, 



(p < 1) 

(33 a) 

J ;1 (p,°) - J\ , 

(p < 1) 

(33b) 

o 

u 

o 

CL 

'w' 

O 

ft 

1 

(p < 1) 

(33c) 


(p < 1) 

(33d) 


A small region which still remains troublesome with the above 
formulas and procedures is that for which £ is small and p is little less 
than unity. In this region, as is shown in the forthcoming paper on evalua- 
tion of the integrals, it is possible to employ a series of the form 



The coefficient integrals of this ^-series surprisingly are evaluable though 
with a little more manipulation than is needed for those of the p-series, so 
this £-series too has has been included as an alternate evaluation procedure 
in FRAC3D. The reader is referred to the forthcoming paper for details. 
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Computations of Crack Stresses in an Infinite Body 

The present analysis for crack stresses in an infinite body began 

by finding expressions for all the stress components and displacements in 

n 1 

terms of a single kind of integral, denoted as 1^ N (p,C), so evaluating 
integrals of that form is the key to stress computations. The most general 
plan for evaluating them is to seek them as entire sets for a given p and £, 
beginning with a base set for a few combinations of indices, and then building 
the rest by recursion. The formulas for integrals of the base set vary some- 
what according to the choice of p and C so several alternative calculations 
must be prepared, but even before this some preparations must be made for 
the numerical quadratures. Thus, several stages of computing are involved 
in evaluating the integrals. 

Gaussian quadrature was chosen for the integrals actually integrated 
numerically. This was advantageous not only as being a well-developed method 
of quadrature but also as providing very efficient approximations for the 
most variable factors in the integrands, that is the Legendre polynomials in 
Equations (18), (19), and (25). To implement this method of quadrature, 
values of the polynomials P^(x) were found for all values of N up to a given 
maximum and for the Gaussian abscissae x for three finenesses of quadrature. 

In order to allow m and k in the Equations (6) to become as large as 10 or 
up to m + 2k = 16 and also to allow for truncation of N as M increases in the 
recursion, values of N as high as 50 were chosen; and, in order to compute 
Legendre polynomials of such high order, double precision arithmetic was used 
here, carrying 30 significant digits in order to counteract the loss of up 
to 15 digits in adding terms of the polynomials. Thus, this most exacting 
part of the calculation was done in advance, before any p or £ was assigned, 
and does not need to be repeated unless the fineness of the quadrature were to 
be changed - which seems quite unnecessary. The finenesses that were selected 
used 32 or 64 or 96 points over the range of integration. After testing showed 
that even the 96-point quadrature ran acceptably fast, it was taken as standard 
for FRAC3D, but results based on other finenesses were employed for accuracy 
checks. The evaluations of P^(x) were multiplied by the customary Gauss 
weights for for the same abscissae so now these products can be used as revised 
weights in the integration and no further evaluations of the P^(x) are needed. 
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, With this background, subroutines were prepared to evaluate the 
alternative bases of integrals to begin the recursion. These subroutines 
are as follows. 


Subroutine 

EVAL 

RHOEQO 

RH0GT1 

RH0LT1 


Range of o and f 

C > 0, p > 0 

C > 0, P = 0 

C - 0 , P > i 

C e p < i 


Equations Employed 
(18) to (24c) 
(25) to (26d) 
(27) to (29d) 
(30) to (33d) 


These subroutines were employed in building another subroutine, called RECUR, 

which uses recursion based on Equations (16) and (17) and builds an entire 

set of integrals for a given p and ? for indices n' = 0,1, 2, and M = 0,1,2,..., 

17, and N = 0,1,2 42 for M = 0, or N = M-l to 43-M for M > 0. RECUR 

includes programming to determine which further subroutine should be called 

for constructing the base set of integrals prior to recursion, including also 

provision for two exceptional regions where EVAL might yield poor results. 

One exceptional region is where 0 < p < 0.7 in which the p-series is to be 

used as a substitute procedure, and the other is where 0.7 ^ p < 1 and 

0 < Q < (l-P)/5 in which the ^-series is to be used. Programming for both the 

p-series and the ^-series is included in RECUR. 

The success of this approach to the evaluation of the integrals 
n 1 

I„ XT (p,£) can be seen in that even with the 96-point quadrature for the bases, 
M,N + 

entire sets of 1000 integrals for a given p and ? are computable on a CDC 

Cyber 73 in little over one second and with good accuracy. By contrast, 

direct integration by the same machine using the original form (1) for the 

integral typically took 40 to 50 seconds per integral and provided less 

accuracy. Thus, the streamlined procedure using RECUR seems to speed the 

4 

integration by a factor up to 10 or more and brings the rate into an accept- 
able range. 
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Stresses Due to Normal and 
Tangential Loads on a Half Space 


Formulation of Effects 
of Surface Loads 


In the crack stress analysis being undertaken, one or more analyses 
for surface-loaded half spaces are to be combined with an analysis for a 
crack in an infinite body in order to represent a lesser body with a crack. 

To make a half-space analysis adaptable for this use, it is important that 
the surface load pattern should be quite flexible. Therefore, an analysis 
is chosen here in which three components of load may be specified arbitrarily 
at many points on the surface of a half space, and it is assumed that surface 
loads at intermediate points can be approximated satisfactorily by linear 
interpolation in either of two perpendicular directions. This representation 
of loads makes them continuous, and this is beneficial in that it avoids the 
artifical singularities of stress that discontinuous loads would imply. 

The points where the load components are specifiable are taken as 
corners of a grid which is rectangular but does not need to be uniform and 
may be incomplete (that is, some boundaries may cross only part of the 
pattern). The loads can be assigned arbitrarily, of course, only at points 
where four rectangles meet since the linear interpolability denies arbitrari- 
ness elsewhere, but this is not a serious limitation. Within this general 
plan, patterns can be drawn which should have ample flexibility to fit a 
wide variety of load distributions. To keep the loads continuous over the 
entire plane, they should of course decrease to zero before they terminate. 

Initially, the elemental load pattern was taken to be the distribu- 
tion of a single component (whether normal or shear) over a rectangle, as 
determined by four arbitrary corner values. This pattern produces stress 
singularities which should be cancelled by singularities caused by combining 
effects of matching loads on adjacent rectangles, but it would be tedious 
to account for this in treating loads on many rectangles. Therefore, an 
alternative kind of elemental load pattern was chosen. This load pattern 
can be related to the earlier load pattern in two stages: first by splitting 

patterns of the earlier kind into four parts of the same general form but 
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with only one non-zero corner value, and second by combining such parts from 
the four rectangles sharing one corner, as illustrated in Figure Al. Elemental 
loads of this kind can be added together to give essentially the same overall 
patterns as would be representable by the earlier rectangular elemental loads, 
but with the newer elements the stress discontinuities are eliminated. These 
new elements are tentatively being called pyramidal since this roughly de- 
scribes their shape, although the faces really are hyperboloidal. 

To formulate stresses induced by the rectangular or pyramidal 
patterns, rectangular coordinates (x,y,z) are introduced with z = 0 being the 
surface of the half space, and auxiliary variables a and 6 are used to represent 
x- and y-distances from a point where stress is to be evaluated to where load 
is applied. In particular, if (x c ,y c> 0) is a corner point of an elemental 
load pattern generating stresses at (x,y,z), then the corresponding corner 
values of a and $ are 

a c “ x c " x » B c " y c " y • (34) 

In the analysis, certain functions of a, 3, and z appear repeatedly, needing 
evaluation for corner values of a and 3, and a certain combination of corner 
values associated with any pyramidal load also appears repeatedly. The 
functions associated with stresses are here denoted as H^(a,3,z) with i iden- 
tifying the load component and j the stress component. In particular j = 

1 (or xx) ^ a , 2 (or yy) 'v. a , 3 (or zz) ^ a , 4 (or yz) ^ t ,5 (or zx) “v 
x y z yz 

T zx* ^ (° r ^ T xy* l° a< ^ component identified by i = 1 is a normal 

load p ^ -a , that for i = 2 is shear load s ^ t , and that for i = 3 is 
z zx 

another shear load t ^ x . (Note that the order of the shear loads differs 

yz 

from that for shear stresses, this being an inheritance from earlier reference 
material.) The functions associated with displacements are denoted as 
yf ^ (pi , (3 , z) , with i used as before, and with j = 1 for u, 2 for v, and 3 for w. 
The functions H 3 (a,3,z) and # 3 (a ,g , z) , to be listed later, were derived as 
shown in Reference 1. The combination of corner evaluations which appears 
repeatedly for any one of these functions (or is the following. 
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sum of four surrounding quarter patterns 



a. Breakdown of rectilinear load pattern into quarter 
patterns determined by the corner values 



*XJ t *xu 


b. Pyramidal load pattern arising from four quarter patterns 
determined by the same corner value 

figure ai. pyramidal load patterns providing load representation 
equivalent to that provided by rectilinear patterns 
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[ H J ■ 73 _h i,iu •(r: + T-h L «i,.u + r-V-»i,uu + 

xX yu 9 xi vll VII > YU VII * 


xi xu yu 


xu yu 


■ 4 (4*4 >m- + (4 + 4X4 + 4^ i * -“ * 4 (4 + 4^ i '“ + <35 > 


+ 1 H . f_l_+ J-V-L- H + 

A x*V i,H WV i,mA 


*£ 


A xuV l * ul 


Here 4^, ^ xu > an< 3 ^y u are x ~ or y - l en 8 t hs of the lower (£) or upper (u) 

rectangles beneath the pyramidal load, and extra subscripts i, m, u applied to 

H^(or call for evaluation using lower, middle, or upper corner values of 

a and ft. The contribution that a pyramidal load with peak value p , s , or 

m m 

t [that is, -cr(x ,y ,0), t (x ,y ,0), or t (x ,y ,0)] makes to stress or 
m mm zx m m yz m m 

displacement components at any point (x,y,z) is readily expressible in terms 
of these combinations of corner evaluations. Indeed, if influence functions 
K 3 (x,y,z) [or K^( x »y» z )] are defined to be 

K|(x,y,z) = - 2^ [Hj(a,e,z)j“ , k{( x , y, z ) - - ^ » z >]* > < 36 > 


then the three elemental pyramidal load patterns on the four-rectangular area 

around (x ,y ,0) produce these stresses and displacements at (x,y,z): 
m m 


2p,u 


j 


o^(x,y,z) 

Eu j (x t y t z) 

d+v) 


P m K l (x,y>2) + s m K 2 (x » y > z) + t m K 3 (x ’ y>z) ’ 
P m Kj(x, y ,z) + s m K^(x,y,z) + t m Kj(x,y,z) , 


(37) 


where the new superscripts by definition imply 


1 

a 


= a 


x’ 




4 

a 


yz 9 


5 

a 


T zx > 


T xy’ 


1 2 3 

U = U, u *= V, u = w 
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To calculate the total stresses and displacements caused at (x,y,z) by the 

three kinds of overall load patterns on the surface of the half space, it is 

necessary only to sum the contributions from all the pyramidal loads. 

In addition to the quantities a and g, each of which is the 

difference of a corner value (such as x or y ) and a stress location 

. , c c 

coordinate (x or y), the and ^ involve also the quantities 


P 





0, = arctan - arctan ^ 
1 a ap 


0 O = arctan - arctan —■ , 

2 B BP * 


0_ •= arctan Q& , 
* z p 


(38) 


and the logarithms ln(p+a), ln(p+g) and ln(p+z). Examination of the functions 
and shows that wherever the argument of one of these logarithms would 
vanish that logarithm is itself multiplied by another vanishing quantity so 
that the product too would vanish. This disposes of all the singularities 
that would have arisen if the elemental load had not been made to vary 
continuously. 

It remains to list the functions H 3 and They are: 


** - 2zp + ffzln(o+p) + 2gzln(g+p) + + (l-2v)|-^ - cyzln(or+p) + 

- S ~^ L - ln(z+p) - ctfe 2 } 

2 2 

H r C 2 ln <«+P> - 2 ^ ln (B + P) + 3 B z e 3 + (l-2v){^ + 


(39a) 


+ 2 * ' ln(a+p) - orzln(z+p) - gz0 2 j- 

apln(crtp) + z 2 ln(g+p) + <*z0 3 + (l-2v)|- gp + oglnCort-p) + 


„xx 

H- - gp 


+ gzln(z+p) - az $ 2 j 
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■ - 2zp + 2azln(orfp) + pzln(3+p) + <*00 3 + - Bzln(£+p) + 

2 2 . 

+ a 5“ ln ( z+ p) ~ 

* yy 2 f 

H 2 ■= ap + z ln(o+p) - or0ln(g+p) + &z0 3 + (l-2v)|- Q 'p + a31n(3+ p ) + 

+ cyzln(z+p) - pzQjJ- 

1 2 2 2 2 
H 3 V " 2 PP " 2 «0 ln (^P> - AL 2L ~ 1 ln(B+p) + 3o-ze 3 + (l-2v)^ ^ + ^-f 2 ” ln(p+p) 

- pzln(z+p) - az0jj 

2 Z 

H 1 - zp + ap0 3 

77 2 

H2 *= -z lti(a + p) - 8 z 0 3 
H 22 *= -z 2 ln(g+p) * <yZ0 3 

H* y ^ 6zln(a+-p) + azln(0+p) - z 2 0 3 + (l- 2 v){^ + pzln(o+p) + ctz ln(p-fp) + 

12 12 

+ a0ln(z+p) + J + p 

2 2 

H* 7 •» Pp - a0ln(cH*p) + z 2 ln(p+p) + az0 3 + (l-2v){- ln(0+p) + 

+ 0zln(z+p) + aze^l 
2 2 

H* 7 ap + z 2 ln(a+p) - apln(0+p) + gze 3 + (1-2 V )|- In(arf-p) + 

+ cr2ln(z+p) + pze 2 } 

H^ 2 «= z 2 ln(p+p) + az0 3 

vz 2 

*= -pzln(a+p) - orzln(0+p) + z 9 3 

H f ** 2zp - 2azln(cH-p) - 0zln(0+p) - Qf00 3 


(39b) 


(39c) 


;} 

(39d) 


(39e) 
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zv 2 

Hj, *= z ln(orfp) + gzQ^ 

*= 2zp - <yzln(a+p) - 2gzin(fl+p) - ( 39 f) 

zx 2 

H3 *= ~pzln(a+p) - orzln(g+p) + z 


^ ln(orfp) - &pzln(p+p) + pz 2 ^ + ^ + 

- ln(orfp) - Qfpzln(p+p) + a ^ a ^ ln(z+p) - — + 

g3e 2 e z \ 

- 6 + T~7 

(40a) 

Hi “ - ^ 2 2z ^ + Q (0 2 " z2 ) 1 - n (Q^p) + " 3z ^ in(p+p) - 2<ypZ0.j + (l-2v)|- + 

+ ^* + ) i n (g+p) + ££°!_zfi_l i n ( 2 + p ) _ opzBjj- 

B(B 2 t . 3? 2 I lnfatp) . °frW) ln(e+p , + £ „ 3 + (1 . 2v) {. 2Bi + 

+ Sfia + g(ft„ -3 . L . ) . ^(o^p) + (dor -3? ) . ln (p +p ) . ^zlnCz+p) - f-fli + 

- £J!2 +fj3 } 

2 6 J 


v 

1 


V 


v 

2 


v 

*3 


&Z£ 

2 


- arpzln(or+p) - Zi SLzZ-l ln(p+p) + QfZ 2 0 3 

. ln(s+|;) . 6l2aia!l ln(2te) 


+ (l-2v)|^g^ + - agzln(a+p) + 

% 2 2 

a ^ org 6 2 a 2 ©31 

" 6 ~ + 2 J 


(40b) 


- ^ Ipfatp) . a^yj . ll>(fH .p) t 4 e 3 + <l-2v){- + 2M 

+ ln(o+p) + °( a ~ 3z ) ln(0+p) - Q0zlti(z-t-p) - a * 9 1 - ~~ + 

- S s t^M + g(B . -- 3z2 ) In (a+p) + p(o' 2 -z 2 )ln(p+ P ) - 2^02^ + <l-2v){- 
+ JL. + fllffi. ~ z ^ ln(or+p) - r ^°’ 2 ~ 6 ^ ln(z+p) - orpz0 2 j- 


+ 



-t 
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, 2 . 2 . 


2 ._ 2 n 


^1 ’"V + "V 1 ' ^ 1 ln(0+p) + (l-2v){^- - -^ + 

■ filfLp-J. ln( Q + p ) - i n (p +p ) + QfBz0 3 | 


,2 2 


#2 * " °2^ + ^ 2 ” ~ Wa+p) + or0zln(0+p) - gz 3 *^ + (l-2v)|^|- + + 


,2 _2 


• ^g ~ ^ ln(o(+p) - agzln(B+p) + ^ ln(z+p) - g ^01 + 

2 

P 3 92 P z l 
’62) 


2 2. 


(40c) 


# 3 * “ ^2^ + o^l^o+p) + ^ 2~ Z ^ l n (B + p) ■ ® z ^03 + (l-2v)^g^ ” ®P z l n (a+p) + 

- ^4^ ln( 6+ p) - ^ ^ ^ } 


Remarks on Calculation of Stresses Due to Surface Loads 


Calculations of stresses due to surface loads involve the same 
quantities p, 0^, 0^, 0^, ln(p+a), ln(p+B) and ln(p+z) many times, so pre- 
calculating them shortens the work. In addition, several other groupings of 
term occur repeatedly in the functions and 7/ 3 , so in the subroutine used 
for calculating surface load effects (SURFAC) these further groupings too are 
precalculated, but more importantly it is recognized that a corner evaluation 
of an H 3 or made for one pyramidal load may be reused for several other 
pyramidal loads sharing that corner. Therefore, before any influence functions 
or k| are calculated by SURFAC for given stress point, (x,y,z) , the functions 
H^(of,p,z) and (cr , (3 , z) are precalculated for all the anticipated corner values 
of a and £ (that is for all df c = x c -x an <3 P c = y c _ y w h ere (x c# y c> 0) ranges over 
all points where rectangles of the surface lattice have corners). In order to 
employ the correct corner evaluations for each pyramidal load, the user assigns 
index numbers to all the corner positions in the plane and then specifies as 
computing input both the position of the pyramidal peak and the indices of the 
nine corners associated with it. The order of specification moves from low x 
to high x, then from high y to low y. By convention, all the pyramidal peak 
positions are numbered first, then the remaining lattice points are numbered. 
(Since it is tedious to assign these index numbers, a computer program LATTICE 
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was prepared to do it.) When the computer receives these specifications for 

the pyramids plus coordinates for the remaining lattice points, the subroutine 

SURFAC finds the required evaluations of and for the given stress point 

and combines them to find the influence functions K/? and K”?. In another mode 

1 1 

of operation, if all the peak loads p^, and t^ are specified, the program 
can compute the stress and displacement components. 

Coordination of Contributing Stress Analyses 

Stress analyses have been formulated for a semi-infinite body 
with an arbitrary set of quasi-pyramidal loads (both normal and shearing) 
acting on its surface, and for an infinite body with a circular crack subject 
to normal and shearing loads expressed by series with arbitrary coefficients. 
Now as an example it will be considered how one may combine three such 
analyses (two surface analyses and one crack analysis) in order to determine 
stresses in a uniformly thick slab with a circular (or part circular) crack, 
with preassigned load (or displacement) distributions on the two faces and 
the crack. A first step in doing this is to determine how coordinates, 
displacements, and stresses expressed in the coordinates for any one of these 
contributing analyses can be reexpressed in terms used in either of the other 
analyses. Then it becomes possible to establish overall boundary conditions 
which imply what hypothetical surface loads should be attributed to each of 
the three loading systems. Evaluation of these hypothetical loads (that is 
evaluation of surface and crack load constants) makes complete stress evalua- 
tion for the slab possible, including stress intensity factors at the crack 
front. 

Coordinate Relationships 

The geometric relationships between the three coordinates systems 
used for the contributing analyses are shown in Figure A2. The rectangular 
system (x,y,z) associated the front face has its axes parallel to those of 
the back-face system (x',y f ,z f ) with the z- and z T -axes colinear but oppo- 
sitely oriented. The x- and x f -axes are oriented alike, but the y- and y f - 


FIGURE A2. COORDINATE SYSTEMS FOR CONTRIBUTING 
STRESS ANALYSES FOR A PLATE 
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axes are not so that both systems can be right-handed. The cylindrical 
coordinates (r,0,p associated with the crack have the plane where 0=0 
coinciding with the xz-plane and the j-axls is directed outward at an in- 
clination m from the positive direction of the x-axis. The crack front has 
r = a and ^ = 0, and the radial offset of its center from the front face of 

the slab is called r , this being positive if that center is outside the 

c 

slab, negative is it is inside. The thickness of the slab is called or T. 

The expressions for back-face coordinates in terms of front-face 
coordinates are 


* 5 

-y > 


Z t - z 


and the inverse expressions are 


x 

y c 

z = 


x 


2 t - z 


(41) 


(41’) 


The front-face coordinates in terms of the cylindrical coordinates can be 
shown to be 


x es r cos 0 sin a) + J cos oo - r c sin u) , 
y = -r sin 0 , 

z = r cos 0 cos a) - ^ sin to - r c cos u) 


(42) 


The expression for the inverse of this relation is chosen to make -rr < 0 < tt, 
and with this provision becomes 


r, ^ .2 2"ll/2 

r «= | (x sin u) + z cos u) + r c ) + y 


0 = arctan 


Zl 

x sinew + z cos u) 


— - 5 < l - E l> E 2 . 


(42') 


^ = x cos u) - z sin (d , 
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where 


= sgn (x s inw + z cos uj + r £ ) 


sgnX ={ X/ [ X I 


-sgn(-y) 

if 

if 


X j- 01 
X 


oj 


and the arctangent takes its principal value (between -rr/2 and tt/2). If the 
denominator of its argument vanishes, the arctangent is put equal to 
tt/2 sgn (numerator) . The relationship between the coordinates for the back 
face and those for the cylindrical system is evidently 


x 7 « r cos 0 sin uj + j cos uj - r c sin uj , 
y 7 ■ r sin 0 , 

z 7 »-r cos 0 cos uj + sinu) + r c cos uj + z t , 
or the inverse form, 


— 2 / 21 
(x 7 sin uj - z 7 cos ou+z t cosu)+r c ) +y I 


1/2 


0 * arctan — — : 


x siruu-z cosuu+z^cosuj+r 

L C 


- \ d-Ej) E 7 , 


where 


j ■ x 7 cos uj + z 7 sin uj - sin uu , 

E 7 ■ sgn (x 7 sinuj - z 7 cos o)+ z fc cos (u+ r c ) , 


(43) 


>(43’) 


E 2 * “ sgn y' , 


with the definition of sgn X as before and with the same principal-value 
choice for the arctangent. If x 7 sin u) - z 7 cos co + z t cos to + r c = y 7 *= 0, 
then 0 is again indeterminate but immaterial since r = 0. 

Stress and displacement components defined in any one system also 
need to be expressed in either of the other systems, so notation is needed 
which identifies the coordinate system as well as the component. Other 
entities often needing identification include the point where the stress is 
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evaluated (the "stress point"), the face where the responsible load element 
is applied, the kind of load component applied, and the load position. To 
provide these identifications for both stresses (a) and displacements (u) , 
superscripts are used to identify in order the stress coordinate systems, the 
stress component, and the stress point, and subscripts are used to identify 
in order the loaded face, the load component, and the load position, but 
unneeded or undetermined indices may be dropped. Entities associated with 
the front face, the back face, or the crack face will be identified respec- 
tively by a, b, or c, and matrix notation will be used to denote sets of 
components. Thus, the values of stress components expressed respectively in 
the front-face, back-face, or crack-face systems are denoted as 






(44) 


and the values of the displacement components are denoted as 



(Here u , u , u are the same as the components u, v, w in the front-face 
x y z 

system: u* , u* , u' are the same as u, v, w in the back-face system; and u , 
x y z a . r 

u Q , U 4 are the same as u, v, w in the crack system.) Evaluation of a * for 

0 * th 

example at a particular stress point, say the n , will be designated by a 

third superscript, making it o a ’^’ n . In similar fashion, a first subscript 

(a, b, c) identifies which face is loaded to produce the stress (provided the 

loading is on a single face). A second subscript shows which kind of load 

is acting. (For loads on rectangular faces, this is 1 for normal loads 
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p 'v -a or -a , , 2 for shearing loads s ^ x or x , , , 3 for shearing loads 
z z zx Z X 

t ^ x or x f , • For loads on the crack face it is 1,2,..., 6 for loads 
yz y z 

associated with the coefficients a 1 , y f , a f ,3 f , y’ , respectively.) A 
third subscript may be used to show which one of a set of load positions 
(on a plane face) or load coefficients (on the crack face) is being used. 

A "load position" on a rectangular face is understood to imply both the peak 
position and the four base lengths used in defining the pyramidal load 
element. The load coefficient designation for the crack accounts for the 
values of both the summation indices m and k. 

Using this new notation, the relationship between stresses and 
displacements in the front-face and back-face systems becomes in matrix 
notation: 

[® b ] ‘ M b,a!y] and 

[»*] * M .,jy > ] and 

1000 o o 
oioo o o 
ooio o o 
0 0 0 1 0 0 
0000-1 0 
0000 0-1 

10 0 
0-10 
,0 0-1 

The relationship between the components in the cylindrical system and those 
in the front-face system can be shown to be as follows. 

M * M c.aM and ’ *0/] • 





[ U J -\,aM 

M - vP] 


(49) 
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where 


H 


. 2 

cos Qsin u) 

2 2 
sin 8 s in uu 

2 

COS U) 


s in ^-0 
2 

cos 0 


- “ sin^sin 2 ua 

cos 0 sin 2 o) 

- sin 20 sin^u) *|sin 20 


0 

0 

0 


2 2 

cos Qcos oj 

2 2 

sin 0 cos tu 

. 2 
8 in a) 


8 in 0 sin 2 aj 

cos 0 sin 2 ai 


-sin 20 cosuj cos 0 sin 2 uj 
2 

sin 20 costu sin 0 sin 2 a) 

0 -sin 2 a) 

cos 0 sinu) -sin 0 cos 2 cju 

sin 0 sina) cos0cos2uj 


-sin 20 sinuu ] 
s in 20 s inu) 

0 

“ COS0COSOJ 
- sin 0 coso) 


12 1 
-- 2 sin 20 cos a) -cos 20 cosu) - — s in 20 s in 2 uo *-cos 20 sinu J 


(50) 


and 




c,a 


cos0sinu) 

-sin0sintij 

COSO) 


-sin0 

“COS 0 

0 


COS 0 CO S(J 0 

-s in0cosa) 
-s in 03 


Inverting these relationships shows 

lyj’VM and M 


(51) 


(49') 


where 


M 


B M 


-1 


a.c c.a 


' 2 , 2 
cos 6s in cu 

2 2 
sin 8s in <ju 

2 

COS U) 

-sin0sin2u; 

cos0s in2a; 

2 

-sin20sin 

sin^Q 

2 

cos 0 

0 

0 

0 

sin20 

2 ft 2 
COS 0COS a) 

2 2 
sin 0cos u) 

. 2 
sin u) 

sin8sin2aj 

-cos0s in2cu 

2 

-sin20cos a) 

- s in20cosaj 

*|sin20cosa) 

0 

cos0s iaa) 

sinQsinou 

-cos20cos(ju 

cos^0sin2u) 

*|s in^0sin2u) 

- s in2u) 

-s in0cos2a) 

cos0cos2u) 

- •- sin20sin2(jo 

- sin20sinuj 

-“sin20sinu) 

0 

- COS0COSU) 

-sin0cosu> 

-cos20sinu) 


and 


n 


a,c 


% 


cos0sino) 
-sin0 
L COS 0 cost!) 


-s in0s into 
-COS 0 
-sin0cosa) 


coscu 

0 

-sino) 


(50') 


(51') 
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It follows that the relationship between 
and those in the back-face system are 

components in the cylindrical system 

[ qC ] ■ M c,bL’ b ] and 

[ u °] ^oPJ • ° r 

i/j ^ciy] and 

(52) 

[ u ] -vM • 

where 


M , = M M 

c,b c,a a,b ’ 

m c,b-\A,b ■ 

M, = M = M M 

d,c c,b b,a a,c 

(53) 

- tff 1 , ** K % 

b,c c,b b,a 'a,c 


Statement of Boundary Conditions 

In order to fit the three contributing stress analyses together 

into a single analysis for a slab with a crack, load coefficients are sought 

for the three analyses so that taken together they yield a least-squares fit 

to boundary conditions prescribed at many points on all three surfaces. To 

state these boundary conditions, it is necessary first to determine the 

influence functions for each of the desired coefficients at all the points 

where the boundary conditions are to be applied. Formulas for the influence 
i i 

functions K^(x,y,z) and K^(x,y,z) for the surface load constants for a single 

load position have been stated already [see Equations (36) to (40)]. The 

same kinds of functions apply, of course, for front and back surface loads. 

Using the revised notation, those influence functions for the front face 

3 - "1 3 1 

loads will now be denoted as K . and K . , since they also express the 

a,i a,i th 

effects in front-face coordinates; and if the load is at the Z position 
and the effects are evaluated at the n ta stress (or boundary) point, then 
they will be called or The influence functions for the back- 

face load constants for the E' 1 "* 1 position, acting on the n^ stress point and 
expressed in back-face coordinates, are given by the same equations but are 
here denoted as or Influence functions for crack load 

constants are readily derived in the crack coordinates from the Equations (6) , 
being products of an F and cos m0 or sin m0 with appropriate sign. Here the 
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influence on the stress component at the n t ^ 1 point due to the l^ 1 kind 

of load coefficient (with i = 1 for a f , ,...,6 for y , ) is denoted as 
. n,k m,k 

c i n 7 7 

*c i*#,"’ with = mn^ + k + 1 where n is the number of k’s used. The 
99 c i n 

influence functions K |t for the displacements are also evident from the 

C 9 1 9 

Equations (6), (Note that yA v’ ^ an< * Yn 1 , will eventually be dropped.) 

0,k 0,k 0,k 

The load constants p , s , t for the £ position on the front 

m m m 

face are now called c c 9 , and c 9 respectively, while those for the 

1 , £ L , Jo 3 , £ 

£ T position on the back face are called c^ ^ t , c’ and c^ • The 

constants ya f . , yg 1 - , ...,yy f . for the crack are now called cV with 
m,k m,k m,k i,£ 

i = 1,2,..., 6 respectively and with £" = mn + k + 1, where n is the total 

tv K 

number of k’s being used. Adding the effects of all front-face loads (which 
include three kinds for each £ = 1,2,...,L), Equation (37) implies that the 

overall components of stress and displacement at the n stress point are 

L 3 


a a ’ j ’ n = Y 7 K a ’j’ a c. # , 

a L L a,i,ji i,A * 


1=1 i=l 
L 3 


a.j.n „ X Y Y v a > J > n c 
u a 2p, L. a,i,i. i,J & ’ 


j = 1,2,, ..,6 , 


j = 1,2,3 


(54) 


(55) 


£* 1 i*=l 


Using matrix notation, Equation (54) may be written as 


_ n n n n 

[ CT a] "[ K a] Icl ■ 


(56) 


where 


n - I n 

0 


r a , l,n 
a 

T a » 2 , n 
a 

j 3 , 3,n 
a 

T a ,4,n 

a 

7 a , 5 , n 
a 

a,6,n 


< c 


> . 


[c] 


1,1 


' 2,1 


'3.1 


' 1,2 


( C 3,L 


(57) 
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and 



K. a 

a, 1,1 

V a > ^ ® 

a, 2,1 

K a ,l,n 
a, 3,1 

i 9 ^ y n 

a, 1,2 

v a 9 ^ 9 ® 

a, 2,2 

- ^ 9 ^ 'I 

a, 3,1, 

g— — i n 

M ■ < 

^a , 2 , n 

a, 1,1 

a , 2 9 n 
K a,2,l 

K a , 2 ,n 
a, 3,1 

, 2 , n 
a, 1,2 

if a , 2 , n 
K a,2,2 

_ _ _ > 2 > a 

a,3,L 1 


K a ,6,n 

a »M 

v a > 6,n 
K a,2,l 

v a > ^ ^ 

a, 3,1 

v a » 6,n 

V,2 

V a 9 ^ 9 n 

K a,2,2 

. t. 3 , 6 , n 

a , 3 , L J 


Equation (55) for the displacements can be written as 



(59) 


where 



a, l,n 
u * * 
a 

J u a ’ 2 ’ n > 
a 

a,3,n 

a 


(60) 


and 



< 


X 

K 

H 


a,l,n 
a, 1,1 

V a 3 ^ 9 ^ 

"a, 2, 1 

v a » 1 > 
a, 3,1 

y a >l,n 

K a,l,2 

y a , l,n 

a, 2, 2 

y a , 1 » ^ 

‘ a, 3,L 

a , 2,n 
a, 1,1 

u a j 2 ? n 
K a,2,l 

y a > 2 > n 

a, 3,1 

v a > 2,n 

a,l,2 

K a,2,n 
“a, 2,2 

y a , 2 , n 
K a,3,L 

a,3,n 
a, 1,1 

v a,3,n 

*a,2,l 

a, 3, n 
X a,3,l 

y a , 3, n 

K a,l,2 

y a »3,n 
X a, 2,2 

- - _ y 3 > 3, n 
a,3,L 


> 


(61) 


In like manner, the overall stresses and displacements caused at the n 
stress point by all the back— face loads (including three for each Z — 1,2,..., 
L') are, if expressed in the back-face coordinates. 
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K] n -KT“'> • 
R]”-i KT'='i • 


(62) 


where |k^ is like the matrix in Equation (58) and is like that in 

Equation (61) except that b's replace all the a's and there are 31/ constants 


'i,£* 


instead of 3L constants c 

th 


i,£ 


Again, the overall stresses and dis- 


placements caused at the n 11 stress point by all the crack loads* (including 
six for each £" = 1,2,...,L") are, if expressed in the crack coordinates, 


[•3 D -MV’ 


(63) 


where is like the matrix in Equation (58) and is like that in 

Equation (61) except that c's replace all the a's and there are 3L" columns 
instead of 3L. Also [ c" ] is like [c] except that it contains 6L" constants 


c! 


i,£" instead of 3L constants c 


i,r 


Of course, the coordinates of the n 


th 


stress point must be evaluated in the system being used, and the Equations 

(41) through (43') make this possible. 

Before the stresses and displacements expressed in different 

systems can be combined, they must be reexpressed into a single system. This 

can be arranged using Equations (46) through (53), but the matrices M , 

c,a 

, and so forth must be evaluated at the chosen stress point, so that they 

C j 3 

are 


denoted as £m c J n , j>y| c ^j n , and so forth for the n 1 "^ stress point. Thus, 
stresses from back 
at the n w “ stress point are** 


the overall stresses from back-face loads expressed in the front-face system 
th 


* If n is the number of m's and n, is the number of k's used for each m there 
m k 


are 6L" crack constants, where L" = n n^, but 3n^ of these are to be dropped 
**0f course J 


= M , , but that kind of simplification does not hold for 
a,b 


matrices involving the crack coordinates. 



A- 39 


p -tti _ — iH _ « — i n — i n _ i — ■ n 

Rl -K, b ] [» b ] -Kb] H Ic,) 


so that 




where [k“] - Kb] [k£] 


(64) 


The overall stresses from crack loads expressed in the front-face system 
similarly are 



where [k*J - Kc] [ K c] 


. (65) 


i-L 

Thus, the overall stresses at the n stress point, being expressed in the 
front-face system are 


»— — j in — i n f— — j n — n 

M '[ K a] [c] + [ K bJ t '' 1 + [ K e] lc " ] 


( 66 ) 


If these stresses are expressed in the back-face system, they become 






(67) 


where 






In the crack- face system, they become 




] + 



» 


( 68 ) 







where 
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In like manner, the displacements at the n stress point, being 
expressed alternatively in the three coordinate systems are 


where 


t uC ]° * i [<f w + i [«b] n ^'i + f M ‘‘" 1 

KHvTKT • KT-kJW 

[ K a] [ M b,al [ X a] • [ X c] * [ M b,c] [ X c] 

KT ■ i\S [tf ■ [ x 'b]" * [»c,br Kf 


(69) 


Equations (66) through (69), together with (58) and (61) and the 

associated definitions of an< ^ so forth, provide formulas for 

computing the overall stress and displacement components if the constants 

an< ^ c i £ " are known. If the constants are not known but values 

of the boundary forces are specified, at least in a least squares sense, at a 

sufficient number of boundary points, then the portions of these equations 

referring to transmissible boundary forces provide means for evaluating those 

constants. It has been arranged that the boundary forces transmissible across 

any face correspond to stress components with j = 3,4,5, when the stresses 

are expressed in the coordinates for that face. (Thus across the front face 

3 . 3 3 4 3 5 

the trensmissible forces correspond to a ’ 'v a , a 5 ^ x , a 9 ^ t ■) 

2 yz 4 a 

If one assembles the influence functions for all front-face loads in producing 
stresses transmissible across that face at all the specified stress points on 
that face (say for n = 1,2,...,N), he has 
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- 

- - 
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R 8 f 4, 2 
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- 

- * 
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K a,5,2 
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R 8 

a , 1 , 2 

R 8 

a, 2, 2 
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- 

- - 
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t.a,5,N ..a,5,N „a, 5,N »a,5,N «a,5,N ,»a,5,N 

«,1, 1 a, 2,1 a, 3, 1 K a,l,2 K a,2,2 K a,3,2 


. K a,5,N 
a,3,L 


(70) 


Similar matrices can be assembled for the influence functions for back-face 
loads acting produce forces transmissible across the front face, that is 

and for the remaining combinations of load face and transmission 
face. Let the stress points on the back face be numbered n' = 1,2,...,N', 

while those on the crack face are numbered n" = 1,2 N". Then the entire 

collection of boundary conditions to be satisfied in the least squares sense 
is 


K: j mI K;m : H k cVm"} 


s 

> > 


s 


[c] 


{a a - j - n } 


[c'] 


{a b ' J ' n '} 


[c"l 


W-"'} 


< / 


V. J 


( 71 ) 


where the specified boundary forces on the three faces are 
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{o 8 ’ j ‘ n } 


r a 8 ’ 3 * 1 ' 


V ’ 3 ’ 1 ' 


^C, 3 , 1 ' 

o 8 * 4 ’ 1 




o C,4>1 

a 8 * 5 * 1 


oW 


a c ’ 5 ’ 1 

a 8 * 3 * 2 


a b > 3 . 2 


a C ’ 3,2 

a 8 * 4 * 2 

> , {c"- 1 -"} • < 

| 


o c,4> 2 > 

a 8 ’ 5 ’ 2 


a b » 5 > 2 


g c ,5,2 

.o 8 > 5 ’ N . 




c,M" 
o > 


(72) 


It is understood here that the indices in each submatrix run through their 

complete range as required, namely three values for j, N for n, N* for n' , 

N" for n", 3 or 6 for i, L for Z, L' for l', and n n. for (Here n is the 

m K m 

number of m's used in Equations (6) and n^ is the number of k’s.) Thus, the 

number of equations implied by Equation (71) is 3N + 3N' +3N", and the number 

of constants to be found is 3L + 3L' + 6n^n^ hut 3n^ of these will be dropped. 

Solving the Equations (71) in the least-squares sense provides 

the constants c, , c*. , , c'.' „ needed to determine stresses and/or displace- 

lyJv 1 j X/ 1 , 36 

ments anywhere in the body using the Equations (66) through (69) or to compute 
crack stress intensity factors using the Equations (11) through (15). 

Computational Procedures 

To implement the above plan for analyzing stresses in a slab with 
a crack, the computing program called FRAC3D has been prepared. What it does, 
briefly, is to accept input data describing the geometry of a slab with a 
crack, together with front-face and back-face lattices, and to call the 
subroutines SURFAC and CRACK for the influence functions related to all the 
assignable load parameters for all the points at which boundary conditions 
are to be imposed. In doing this, FRAC3D calls the subroutine TRANSF to find 
the needed coordinate transformations using Equations (41) through (43'), and 
it calls the subroutine CONVRT to find the needed stress transformations 
by using Equations (46) through (53). It coordinates the additive 
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stresses as in Equations (67) through (70) so that contributions from all 
three loading systems are expressed in the coordinate system to which the 
particular boundary conditions refer. Then it assembles all the results into 
a matrix of the form displayed in Equations (70) through (72), readying them 
for solution for the loading constants using a previously available program 
called MATSOL. 

The program FRAC3D is capable of treating geometries and loadings 
having wide generality, and in particular, having little or no symmetry. If, 
however, the crack is perpendicular to the faces of the slab and the loading 
has symmetry around both the xz- and yz-planes, then the matrix is unduly 
large. For such cases, a subroutine AXYZ has been included which combines 
columns of the matrix which relate to unknown constants which must be equal 
or merely of opposite sign on account of the symmetry. As a result of this 
reduction of columns, the rows of the matrix for stress points symmetrically 
situated with respect to the coordinate axes on each face become identical. 
Therefore, boundary conditions for front and back faces are assigned for only 
one quadrant (including its edges), and those for the crack only for 0 >_ 0. 

The program weights each row according to the number of boundary conditions it 
represents, the weights being 1,^/2 or 2 according to whether there are 1, 2 or 
4 represented. This yields a much reduced system of equations to submit to 
MATSOL, but the computed results are the same as if the matrix has not been 
reduced. This capacity of AXYZ actually extends to allowing for several 
kinds of symmetry and to treating symmetry when the body includes as many as 
six sides of a rectangular parallelopiped. 

The program MATSOL performs a least-square solution for the 

loading constants for both flat faces and the crack. Once these have been 

found, in particular once the crack constants a' , , S’ ,,..., are determined, 

m , k 

it is straightforward to find the variation of stress intensity factors, 
along the crack, using Equations (11) through (15) or to get half the crack 
opening displacement, using the surface displacement w as computed by CRACK 
from one of the Equations (6) . The program also computes residuals at the 
points where boundary conditions are fitted and provides constants which 
could be applied in the Equations (6), (39), and (40) to find stresses or 
displacements anywhere in the slab. The program MATSOL has an option also 
of using several exact equations in conjunction with those being solved in 



A-44 


the least squares sense. Such an exact equation can be used to impose an 
equilibrium equation among the applied loads, and that is usually done as a 
matter of propriety. An expanded version of MATSOL has editing capability 
by which blocks of equations may be weighted differently and selected load 
constants can be removed from being used. (It is not necessary to use MATSOL 
to remove the spurious crack constants Yq ^ and Yq however, since 

FRAC3D does that automatically.) 

Use of the program FRAC3D requires choosing of surface lattices, 

choosing the structure of crack function series, and assignment of boundary 

conditions. Such matters are discussed in the body of this report and also 

in the report to which this manual is a supplement [2], It may be noted here 

that some structuring of the crack function series can be accomplished by 

using some simple, available inserts for the subroutine AXYZ which can remove 

terms for odd m's and then also terms for which m4-2k > 2n. , if that is 

k 

desired and if the only constants are The early structuring is handled 

through the assignment of the n^ and n^ to be applied for each of the six 
kinds of crack load constants. 
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APPENDIX B 

COMPONENTS OF FRAC3D 


The program FRAC3D has many subroutines which are employed according 
to the location of the points and kinds of stresses being considered at any 
stage of the calculation. A brief description of the various parts follows. 

Program FRAC3D 

This is the part of the program which reads input and starts a course 
of action depending on the mode, on the face where a boundary condition is 
assigned, and on the system in which the load under consideration acts. 

Detailed calculations are referred to appropriate subroutines but returned 
to FRAC3D for reporting as output. 

In the SOLVE mode, the intent of FRAC3D is to construct a matrix of the 
form shown in Equation (71) of Appendix A. The notation in that equation, 
and discussion of how its parts are formed are provided in the subsection of 
Appendix A containing that equation. In the RESULT mode, similar matrix rows 
are formed for each stress evaluation point, with all stress components being 
retained and all displacement components being included. The elements of 
each row are multiplied by the corresponding load constants and summed to 
determine the stresses and displacements. 

Brief descriptions of the subroutines follow: 

SURFAC - This calculates influence functions and other entities related 
to or arising from loads on the plane surfaces. Equations used 
by this subroutine are shown as Equations (35) through (40) in 
Appendix A. It is presumed the input specifies lattice points 
and pyramidal base indices as provided by the program LATTICE. 

CRACK - This directs calculations for influence functions and other 

entities related to or arising from loads on the crack. In its 
assembly work it employs Equations (3) through (9) of Appendix A. 
In doing this, it refers evaluation of the many needed integrals 
to other subroutines via the subroutine RECUR. 
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AUXEQ - This assembles the auxiliary equations relating to overall 

equilibrium. These equations are employed in the SOLVE mode, but 
are kept apart from the regular boundary conditions. After form- 
ation of normal equations from the boundary-condition equations 
by MATSOL, these equilibrium equations are attached to the system 
by a bordering process. Thereafter, they are simply members of 
the system which MATSOL solves to find the load constants. 

CONVRT - This performs conversions of stress components from one coordi- 
nate system to another. Equations for the conversions are given 
by Equations (44) to (53) of Appendix A. 

TRANSF - This performs transformations of coordinates between cylindrical 
and rectangular coordinate systems. The equations used for this 
purpose are provided by Equations (41) to (43* ) of Appendix A. 

AXYZ - This constructs the symmetry chain which is described on Pages 
14 to 16 and illustrated by Pages 50 to 52. This routine on a 
logical, step-by-step process with many branches. 

RECUR - This is the controlling subroutine for evaluation of the inte- 
grals used by the crack functions. Usually it assigns some other 
subroutine the task of evaluating a base set of the integrals, as 
described in Appendix A, and on their return it applies recursion 
processes to enlarge the set of integrals greatly. The recursion 
is two dimensional, being based on Equations (16) and (17) of 
Appendix A. Two special classes of integrals are treated almost 
entirely within RECUR itself, namely those for integrals with 
p < 0.7 (treated by the p-series) or with both 0.7 $ p < 1 and 
0 < £ < (l-p)/5 (treated by the ^-series). These special series 
are only sketched in Appendix A, but they derived and discussed 
rather fully in the forthcoming paper mentioned in connection 
with those series in Appendix A. The use of recursion following 
efficient evaluation of base sets of integrals explains the great 
speed of evaluation of the crack integrals, but recursion can 
magnify numerical errors rapidly. As the forthcoming paper ex- 
plains, accuracy as well as speed was studied until a good pro- 
cedure was found for RECUR to select for any combination (p,£). 
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EVAL 

RHOEQO 

RHOGT1 

RHOLT1 

LOGF1 

TH1F 

TH3F 

BLKDAT 


This computes the most frequently used type of base integrals by 
a modified Gaussian quadrature of integrals greatly transformed 
by analytic processes. It employs mainly the Equations (18) to 
(24) of Appendix A. The quadrature is greatly speeded by 
precalculations of many Legendre polynomials multiplied by Gauss 
weights, with the results being recorded in the subroutine BLKDAT. 
This computes base sets of crack integrals for cases with p = 0. 
It employs Equations (25) and (26) of Appendix A. 

This computes sets of crack integrals for cases with Q = 0 and 
p > 1. It uses Gaussian quadrature applied to Equations (27) to 
(29) of Appendix A. 

This computes base sets of crack integrals for cases with Q - 0 
and p < 1 transformed into polynomials as shown by Equations (30) 
to (33) of Appendix A. 

This is a tool used by SURFAC to eliminate indeterminacies in the 
surface load functions. 

This computes a special function arising in SURFAC, denoted as 0^ 
in the Equations (38) of Appendix A. 

This computes a special function arising in SURFAC, denoted as 0^ 
in the Equations (38) of Appendix A. 

This contains data for the Gaussian quadrature performed by RH0GT1 
and for the modified Gaussian quadrature performed by EVAL with 
Gauss weights multiplied by Legendre polynomials. 


